In this project we investigated the electrophysiological correlates of the combined processing of basic visual properties (i.e., size and contrast) and the emotional content of simple words.

GRAPHS

grand.avg.localizer <-
  read_csv(here::here("EEG/main/data/main_grandEEG_localizer.csv")) %>%
  select(-condition) %>% # drop condition column
  gather( # convert to long format
    key = timepoint, # name of new column
    value = amplitude, # value of new column
    -c(participant, electrode) # keep participants and electrodes as columns
  )

# summarized data from each time point (& within-subject 95% CI)
grand.avg.localizer.pointsummary <-
  grand.avg.localizer %>%
  summarySEwithin(
    data = .,
    measurevar = "amplitude",
    withinvars = c("electrode", "timepoint"),
    idvar = "participant"
  ) %>%
  mutate(
    timepoint = as.numeric(levels(timepoint))[timepoint] # re-convert time points to numeric
  )

# plot
ggplot(
  grand.avg.localizer.pointsummary,
  aes(
    x = timepoint,
    y = amplitude,
    group = electrode
  )
) +
  geom_vline( # vertical reference line
    xintercept = 0,
    linetype = "dashed",
    color = "black",
    size = 1.2,
    alpha = .8
  ) +
  geom_hline( # horizontal reference line
    yintercept = 0,
    linetype = "dashed",
    color = "black",
    size = 1.2,
    alpha = .8
  ) +
  geom_vline( # vertical reference lines
    xintercept = seq(-200, 1000, 100),
    linetype = "dotted",
    color = "#999999",
    size = .8,
    alpha = .5
  ) +
  geom_hline( # horizontal reference lines
    yintercept = seq(-3, 3, 1),
    linetype = "dotted",
    color = "#999999",
    size = .8,
    alpha = .5
  ) +
  geom_line( # one line per electrode
    size = 1.2,
    color = "#494847",
    alpha = .8
  ) +
  geom_ribbon( # 95% CI
    aes(
      ymin = amplitude - ci,
      ymax = amplitude + ci
    ),
    linetype = "dotted",
    size = .1,
    alpha = .1,
    show.legend = FALSE
  ) +
  labs(
    title = "", # title & axes labels
    x = "time (ms)",
    y = expression(paste("amplitude (", mu, "V)"))
  ) +
  scale_x_continuous(breaks = seq(-200, 1000, 100)) + # x-axis: tick marks
  scale_y_reverse(
    breaks = seq(-3, 3, 1), # y-axis: tick marks
    limits = c(3, -3)
  ) +
  theme_classic(base_size = 20) +
  theme(
    plot.title = element_text(
      size = 28,
      hjust = .5,
      face = "bold"
    ),
    legend.position = "none"
  )
Butterfly plot

Butterfly plot

Prepare data for:

  • topographies
  • waveforms
  • (RDI plots) pirateplots
grand.avg.comps <-
  read_csv(here::here("EEG/main/data/main_avgERP.csv")) %>%
  gather( # convert to long format
    key = timepoint,
    value = amplitude,
    -c(participant, elec_cluster, condition)
  ) %>%
  # rename levels
  mutate(
    condition = revalue(
      as.factor(condition),
      c(
        "all" = "localizer",
        "negative" = "negative",
        "negative minus neutral" = "negative minus neutral",
        "negativeLargeBright" = "negative large bright",
        "negativeLargeDark" = "negative large dark",
        "negativeSmallBright" = "negative small bright",
        "negativeSmallDark" = "negative small dark",
        "neutral" = "neutral",
        "neutralLargeBright" = "neutral large bright",
        "neutralLargeDark" = "neutral large dark",
        "neutralSmallBright" = "neutral small bright",
        "neutralSmallDark" = "neutral small dark"
      )
    ),
    # variable with main effect of size
    size = revalue(
      condition,
      c(
        "localizer" = NA,
        "negative" = NA,
        "negative minus neutral" = NA,
        "negative large bright" = "large",
        "negative large dark" = "large",
        "negative small bright" = "small",
        "negative small dark" = "small",
        "neutral" = NA,
        "neutral large bright" = "large",
        "neutral large dark" = "large",
        "neutral small bright" = "small",
        "neutral small dark" = "small"
      )
    ),
    # variable with main effect of contrast
    cont = revalue(
      condition,
      c(
        "localizer" = NA,
        "negative" = NA,
        "negative minus neutral" = NA,
        "negative large bright" = "bright",
        "negative large dark" = "dark",
        "negative small bright" = "bright",
        "negative small dark" = "dark",
        "neutral" = NA,
        "neutral large bright" = "bright",
        "neutral large dark" = "dark",
        "neutral small bright" = "bright",
        "neutral small dark" = "dark"
      )
    ),
    # variable with main effect of emotion
    emo = revalue(
      condition,
      c(
        "localizer" = NA,
        "negative" = NA,
        "negative minus neutral" = NA,
        "negative large bright" = "negative",
        "negative large dark" = "negative",
        "negative small bright" = "negative",
        "negative small dark" = "negative",
        "neutral" = NA,
        "neutral large bright" = "neutral",
        "neutral large dark" = "neutral",
        "neutral small bright" = "neutral",
        "neutral small dark" = "neutral"
      )
    )
  )

P1

Amplitude extracted from the following electrode cluster: P7 P9 PO7 O1 O2 PO8 P8 P10.

# summarized data from each time point, including within-subject 95% CIs
P1.grand.avg.comps.pointsummary <-
  grand.avg.comps %>% # data frame
  filter( # filter out unused conditions
    elec_cluster == "P1",
    condition != "negative",
    condition != "negative minus neutral",
    condition != "neutral"
  ) %>%
  summarySEwithin(
    data = .,
    measurevar = "amplitude",
    withinvars = c("size", "cont", "emo", "condition", "timepoint"),
    idvar = "participant"
  ) %>%
  unite(visualfeats, c("size", "cont"), sep = " ") %>% # for graph
  mutate(
    visualfeats = as.factor(ifelse(visualfeats == "NA NA", "localizer", visualfeats)),
    timepoint = as.numeric(levels(timepoint))[timepoint]
  )

# plot waveforms
ggplot(
  data = filter(
    P1.grand.avg.comps.pointsummary,
    condition != "localizer"
  ),
  aes(
    x = timepoint,
    y = amplitude,
    linetype = emo,
    color = visualfeats
  )
) +
  geom_vline( # vertical reference line
    xintercept = 0,
    linetype = "dashed",
    color = "black",
    size = 1.2,
    alpha = .8
  ) +
  geom_hline( # horizontal reference line
    yintercept = 0,
    linetype = "dashed",
    color = "black",
    size = 1.2,
    alpha = .8
  ) +
  geom_vline( # vertical reference lines
    xintercept = seq(-200, 1000, 50),
    linetype = "dotted",
    color = "#999999",
    size = .8,
    alpha = .5
  ) +
  geom_hline( # horizontal reference lines
    yintercept = seq(-3, 3, 1),
    linetype = "dotted",
    color = "#999999",
    size = .8,
    alpha = .5
  ) +
  geom_line( # one line per condition
    size = 2.5,
    alpha = .6
  ) +
  geom_line( # thick black line (localizer)
    data = filter(
      P1.grand.avg.comps.pointsummary,
      visualfeats == "localizer"
    ),
    aes(
      x = timepoint,
      y = amplitude
    ),
    size = 5,
    show.legend = FALSE,
    inherit.aes = FALSE
  ) +
  geom_ribbon( # localizer 95% CI
    data = filter(
      P1.grand.avg.comps.pointsummary,
      visualfeats == "localizer"
    ),
    aes(
      ymin = amplitude - ci,
      ymax = amplitude + ci
    ),
    linetype = "dotted",
    size = .1,
    alpha = .3,
    show.legend = FALSE
  ) +
  scale_color_manual(
    breaks = c(
      "localizer", "large bright", # re-order conditions
      "large dark", "small bright", "small dark"
    ),
    values = waveforms.palette
  ) + # line colors
  scale_fill_manual(values = waveforms.palette) + # ribbon colors
  labs(
    title = "P1 (66-148 ms)", # title & axes labels
    x = "time (ms)",
    y = expression(paste("amplitude (", mu, "V)"))
  ) +
  scale_x_continuous(
    breaks = seq(-200, 1000, 50), # x-axis: tick marks
    limits = c(-100, 250)
  ) +
  scale_y_reverse(
    breaks = seq(-3, 3, 1), # y-axis: tick marks
    limits = c(3, -3)
  ) +
  annotate("rect", # highlight time window used for analysis
    xmin = 66,
    xmax = 148,
    ymin = -.5,
    ymax = 3,
    linetype = "solid",
    size = 2,
    color = "#de1d1d",
    alpha = 0
  ) +
  theme_classic(base_size = 20) +
  theme(
    plot.title = element_text(
      size = 28,
      hjust = .5,
      face = "bold"
    ),
    legend.position = "none"
  )

N1

Amplitude extracted from the following electrode cluster: TP7 P7 P9 TP8 P8 P10.

N1.grand.avg.comps.pointsummary <-
  grand.avg.comps %>%
  filter(
    elec_cluster == "N1",
    condition != "negative",
    condition != "negative minus neutral",
    condition != "neutral"
  ) %>%
  summarySEwithin(
    data = .,
    measurevar = "amplitude",
    withinvars = c("size", "cont", "emo", "condition", "timepoint"),
    idvar = "participant"
  ) %>%
  unite(visualfeats, c("size", "cont"), sep = " ") %>%
  mutate(
    visualfeats = as.factor(ifelse(visualfeats == "NA NA", "localizer", visualfeats)),
    timepoint = as.numeric(levels(timepoint))[timepoint]
  )

# plot waveforms
ggplot(
  data = filter(
    N1.grand.avg.comps.pointsummary,
    condition != "localizer"
  ),
  aes(
    x = timepoint,
    y = amplitude,
    linetype = emo,
    color = visualfeats
  )
) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed",
    color = "black",
    size = 1.2,
    alpha = .8
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed",
    color = "black",
    size = 1.2,
    alpha = .8
  ) +
  geom_vline(
    xintercept = seq(-200, 1000, 50),
    linetype = "dotted",
    color = "#999999",
    size = .8,
    alpha = .5
  ) +
  geom_hline(
    yintercept = seq(-3, 3, 1),
    linetype = "dotted",
    color = "#999999",
    size = .8,
    alpha = .5
  ) +
  geom_line(
    size = 2.5,
    alpha = .6
  ) +
  geom_line(
    data = filter(
      N1.grand.avg.comps.pointsummary,
      visualfeats == "localizer"
    ),
    aes(
      x = timepoint,
      y = amplitude
    ),
    size = 5,
    show.legend = FALSE,
    inherit.aes = FALSE
  ) +
  geom_ribbon(
    data = filter(
      N1.grand.avg.comps.pointsummary,
      visualfeats == "localizer"
    ),
    aes(
      ymin = amplitude - ci,
      ymax = amplitude + ci
    ),
    linetype = "solid",
    size = .1,
    alpha = .3,
    show.legend = FALSE
  ) +
  scale_color_manual(
    breaks = c(
      "localizer", "large bright",
      "large dark", "small bright", "small dark"
    ),
    values = waveforms.palette
  ) +
  scale_fill_manual(values = waveforms.palette) +
  labs(
    title = "N1 (150-260 ms)",
    x = "time (ms)",
    y = expression(paste("amplitude (", mu, "V)"))
  ) +
  scale_x_continuous(
    breaks = seq(-200, 1000, 50),
    limits = c(-100, 400)
  ) +
  scale_y_reverse(
    breaks = seq(-3, 3, 1),
    limits = c(3, -3)
  ) +
  annotate("rect",
    xmin = 150,
    xmax = 260,
    ymin = -2.5,
    ymax = 1,
    linetype = "solid",
    size = 2,
    color = "#de1d1d",
    alpha = 0
  ) +
  theme_classic(base_size = 20) +
  theme(
    plot.title = element_text(
      size = 28,
      hjust = .5,
      face = "bold"
    ),
    legend.position = "none"
  )

EPN

Because the topography differs from the classical EPN topography published in the literature, we identified this component by calculating the difference between negative and neutral conditions (averaged across size and contrast). This is in accordance with our pre-registered protocol.

Note that this was not necessary when identifying the LPP.

Amplitude extracted from the following electrode cluster: P7 P9 PO7 PO3 O1 Oz Iz O2 PO4 PO8 P8 P10.

EPN.grand.avg.comps.pointsummary <-
  grand.avg.comps %>%
  filter(
    elec_cluster == "EPN",
    condition == "negative minus neutral" |
      condition == "negative large bright" |
      condition == "negative large dark" |
      condition == "negative small bright" |
      condition == "negative small dark" |
      condition == "neutral large bright" |
      condition == "neutral large dark" |
      condition == "neutral small bright" |
      condition == "neutral small dark"
  ) %>%
  summarySEwithin(
    data = .,
    measurevar = "amplitude",
    withinvars = c("size", "cont", "emo", "condition", "timepoint"),
    idvar = "participant"
  ) %>%
  unite(visualfeats, c("size", "cont"), sep = " ") %>% # for graph
  mutate(
    visualfeats = as.factor(ifelse(visualfeats == "NA NA", "negative minus neutral", visualfeats)),
    timepoint = as.numeric(levels(timepoint))[timepoint]
  )

# plot waveforms
ggplot(
  data = EPN.grand.avg.comps.pointsummary,
  aes(
    x = timepoint,
    y = amplitude,
    linetype = emo,
    color = visualfeats
  )
) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed",
    color = "black",
    size = 1.2,
    alpha = .8
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed",
    color = "black",
    size = 1.2,
    alpha = .8
  ) +
  geom_vline(
    xintercept = seq(-200, 1000, 50),
    linetype = "dotted",
    color = "#999999",
    size = .8,
    alpha = .5
  ) +
  geom_hline(
    yintercept = seq(-3, 3, 1),
    linetype = "dotted",
    color = "#999999",
    size = .8,
    alpha = .5
  ) +
  geom_line(
    size = 2.5,
    alpha = .6
  ) +
  geom_line(
    data = filter(
      EPN.grand.avg.comps.pointsummary,
      visualfeats == "negative minus neutral"
    ),
    aes(
      x = timepoint,
      y = amplitude
    ),
    size = 5,
    show.legend = FALSE,
    inherit.aes = FALSE
  ) +
  geom_ribbon(
    data = filter(
      EPN.grand.avg.comps.pointsummary,
      visualfeats == "negative minus neutral"
    ),
    aes(
      ymin = amplitude - ci,
      ymax = amplitude + ci
    ),
    linetype = "dotted",
    size = .1,
    alpha = .3,
    show.legend = FALSE
  ) +
  scale_color_manual(
    breaks = c(
      "negative minus neutral", "large bright",
      "large dark", "small bright", "small dark"
    ),
    values = waveforms.palette
  ) +
  scale_fill_manual(values = waveforms.palette) +
  labs(
    title = "EPN (300-500 ms)",
    x = "time (ms)",
    y = expression(paste("amplitude (", mu, "V)"))
  ) +
  scale_x_continuous(
    breaks = seq(-200, 1000, 50),
    limits = c(-100, 600)
  ) +
  scale_y_reverse(
    breaks = seq(-2, 3, 1),
    limits = c(3, -2)
  ) +
  annotate("rect",
    xmin = 300,
    xmax = 500,
    ymin = -1,
    ymax = 2.5,
    linetype = "solid",
    size = 2,
    color = "#de1d1d",
    alpha = 0
  ) +
  theme_classic(base_size = 20) +
  theme(
    plot.title = element_text(
      size = 28,
      hjust = .5,
      face = "bold"
    ),
    legend.position = "none"
  )

LPP

Amplitude extracted from the following electrode cluster: P1 Pz P2 P4 P6 P8 P10 POz PO4 PO8.

LPP.grand.avg.comps.pointsummary <- 
  grand.avg.comps %>%
  filter(
    elec_cluster == "LPP",
    condition == "localizer" |
      condition == "negative large bright" |
      condition == "negative large dark" |
      condition == "negative small bright" |
      condition == "negative small dark" |
      condition == "neutral large bright" |
      condition == "neutral large dark" |
      condition == "neutral small bright" |
      condition == "neutral small dark"
  ) %>%
  summarySEwithin(
    data = .,
    measurevar = "amplitude",
    withinvars = c("size", "cont", "emo", "condition", "timepoint"),
    idvar = "participant"
  ) %>%
  unite(visualfeats, c("size", "cont"), sep = " ") %>% # for graph
  mutate(
    visualfeats = as.factor(ifelse(visualfeats == "NA NA", "localizer", visualfeats)),
    timepoint = as.numeric(levels(timepoint))[timepoint]
  )

# plot waveforms
ggplot(
  data = LPP.grand.avg.comps.pointsummary,
  aes(
    x = timepoint,
    y = amplitude,
    linetype = emo,
    color = visualfeats
  )
) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed",
    color = "black",
    size = 1.2,
    alpha = .8
  ) +
  geom_hline(
    yintercept = 0,
    linetype = "dashed",
    color = "black",
    size = 1.2,
    alpha = .8
  ) +
  geom_vline(
    xintercept = seq(-200, 1000, 100),
    linetype = "dotted",
    color = "#999999",
    size = .8,
    alpha = .5
  ) +
  geom_hline(
    yintercept = seq(-3, 3, 1),
    linetype = "dotted",
    color = "#999999",
    size = .8,
    alpha = .5
  ) +
  geom_line(
    size = 2.5,
    alpha = .6
  ) +
  geom_line(
    data = filter(
      LPP.grand.avg.comps.pointsummary,
      visualfeats == "localizer"
    ),
    aes(
      x = timepoint,
      y = amplitude
    ),
    size = 5,
    show.legend = FALSE,
    inherit.aes = FALSE
  ) +
  geom_ribbon(
    data = filter(
      LPP.grand.avg.comps.pointsummary,
      visualfeats == "localizer"
    ),
    aes(
      ymin = amplitude - ci,
      ymax = amplitude + ci
    ),
    linetype = "dotted",
    size = .1,
    alpha = .3,
    show.legend = FALSE
  ) +
  scale_color_manual(
    breaks = c(
      "localizer", "large bright",
      "large dark", "small bright", "small dark"
    ),
    values = waveforms.palette
  ) +
  scale_fill_manual(values = waveforms.palette) +
  labs(
    title = "LPP (402-684 ms)",
    x = "time (ms)",
    y = expression(paste("amplitude (", mu, "V)"))
  ) +
  scale_x_continuous(
    breaks = seq(-200, 1000, 100),
    limits = c(-100, 800)
  ) +
  scale_y_reverse(
    breaks = seq(-2, 3, 1),
    limits = c(3, -2)
  ) +
  annotate("rect",
    xmin = 402,
    xmax = 684,
    ymin = -.5,
    ymax = 2.5,
    linetype = "solid",
    size = 2,
    color = "#de1d1d",
    alpha = 0
  ) +
  theme_classic(base_size = 20) +
  theme(
    plot.title = element_text(
      size = 28,
      hjust = .5,
      face = "bold"
    ),
    legend.position = "none"
  )



CONFIRMATORY ANALYSIS

Calculate and compare the Bayes Factor of different linear mixed-effects models. Participants’ ID is included as random factor and its variance set as nuisance.

Compare (against the null model) the following models:

  1. main effects of size and emotion
  2. interactive effects of size and emotion
  3. main effects of contrast and emotion
  4. interactive effects of contrast and emotion
  5. main effects of size, contrast, and emotion
  6. interactive effects of size, contrast, and emotion

Afterwards, compare the best competing models to understand which one should be preferred overall.

P1

size cont emo N P1 sd se ci
large dark negative 2374 1.00 3.82 0.08 0.15
large dark neutral 2369 1.06 3.76 0.08 0.15
large bright negative 2373 1.21 3.76 0.08 0.15
large bright neutral 2364 1.34 3.76 0.08 0.15
small dark negative 2383 1.09 3.81 0.08 0.15
small dark neutral 2376 0.92 3.73 0.08 0.15
small bright negative 2365 0.34 3.83 0.08 0.15
small bright neutral 2363 0.38 3.86 0.08 0.16
for (k in 1:length(scaling.factor)) { # loop through scaling factors

  ### main effects of size and emotion
  P1.sizeplusemo.BF <- lmBF(
    P1 ~ size + emo + participant, # formula
    data.EEG.trial, # data
    whichRandom = "participant", # random effects
    rscaleFixed = scaling.factor[k], # scaling factor of prior on effect size
    rscaleRandom = "nuisance", # prior scale for standardized random effects
    rscaleCont = "medium", # prior scale for standardized slopes
    method = "simple", # MCMC sampling
    iterations = niter # number of MCMC iterations
  )

  saveRDS(P1.sizeplusemo.BF, # save model
    file = here::here(paste0("EEG/main/models/P1_sizeplusemo_BF_", scaling.factor[k], ".rds"))
  )

  ### interactive effects of size and emotion
  P1.sizebyemo.BF <- lmBF(
    P1 ~ size * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.sizebyemo.BF,
    file = here::here(paste0("EEG/main/models/P1_sizebyemo_BF_", scaling.factor[k], ".rds"))
  )

  ### main effects of contrast and emotion
  P1.contplusemo.BF <- lmBF(
    P1 ~ cont + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.contplusemo.BF,
    file = here::here(paste0("EEG/main/models/P1_contplusemo_BF_", scaling.factor[k], ".rds"))
  )

  ### interactive effects of contrast and emotion
  P1.contbyemo.BF <- lmBF(
    P1 ~ cont * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.contbyemo.BF,
    file = here::here(paste0("EEG/main/models/P1_contbyemo_BF_", scaling.factor[k], ".rds"))
  )

  ### main effects of size, contrast, and emotion
  P1.sizepluscontplusemo.BF <- lmBF(
    P1 ~ size + cont + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.sizepluscontplusemo.BF,
    file = here::here(paste0("EEG/main/models/P1_sizepluscontplusemo_BF_", scaling.factor[k], ".rds"))
  )

  ### interactive effects of size, contrast, and emotion
  P1.sizebycontbyemo.BF <- lmBF(
    P1 ~ size * cont * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.sizebycontbyemo.BF,
    file = here::here(paste0("EEG/main/models/P1_sizebycontbyemo_BF_", scaling.factor[k], ".rds"))
  )
}
# preallocate matrix with all BFs
compare.P1.BF <- matrix(NA, 6, length(scaling.factor))

# preallocate matrix with all % errors
compare.P1.perc.err <- matrix(NA, 6, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) { # loop through scaling factors

  ### load models
  P1.sizeplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_sizeplusemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.sizebyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_sizebyemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.contplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_contplusemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.contbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_contbyemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.sizepluscontplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_sizepluscontplusemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.sizebycontbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_sizebycontbyemo_BF_", scaling.factor[k], ".rds"))
  )

  # BFs
  compare.P1.BF[, k] <- c(
    P1.sizeplusemo.BF@bayesFactor$bf,
    P1.sizebyemo.BF@bayesFactor$bf,
    P1.contplusemo.BF@bayesFactor$bf,
    P1.contbyemo.BF@bayesFactor$bf,
    P1.sizepluscontplusemo.BF@bayesFactor$bf,
    P1.sizebycontbyemo.BF@bayesFactor$bf
  )

  # percentage of error
  compare.P1.perc.err[, k] <- c(
    P1.sizeplusemo.BF@bayesFactor$error * 100,
    P1.sizebyemo.BF@bayesFactor$error * 100,
    P1.contplusemo.BF@bayesFactor$error * 100,
    P1.contbyemo.BF@bayesFactor$error * 100,
    P1.sizepluscontplusemo.BF@bayesFactor$error * 100,
    P1.sizebycontbyemo.BF@bayesFactor$error * 100
  )
}

# summary confirmatory analysis
compare.P1 <- data.frame(
  "model" = c("size + emo", "size x emo", "contr + emo", "cont x emo", "size + cont + emo", "size x cont x emo"),
  "nar" = compare.P1.BF[, 1], "nar.p.err" = compare.P1.perc.err[, 1],
  "med" = compare.P1.BF[, 2], "med.p.err" = compare.P1.perc.err[, 2],
  "wid" = compare.P1.BF[, 3], "wid.p.err" = compare.P1.perc.err[, 3]
)

# sort according to medium scaling factor (in descending order)
compare.P1 <- compare.P1[order(compare.P1$med, decreasing = TRUE), ]

# nicer table
kable(compare.P1)
model nar nar.p.err med med.p.err wid wid.p.err
6 size x cont x emo 516.892 16.2436 514.126 18.3575 512.777 39.6341
5 size + cont + emo 491.855 13.1956 490.608 10.9562 489.588 13.3064
1 size + emo 487.969 11.3261 487.267 10.8352 486.329 10.6004
2 size x emo 485.359 13.1905 484.269 12.3227 483.004 12.4829
3 contr + emo 455.029 11.1600 454.463 10.1229 453.767 10.9300
4 cont x emo 452.174 20.2480 451.206 14.2266 449.945 11.2667

When using a JZS prior with scaling factor r = 0.707 placed on standardized effect sizes, the model size x cont x emo ought to be preferred.
The best model (size x cont x emo) explains the observed data 1.6369e+10 times better than the second best model (size + cont + emo).

Paired comparisons

# interactions of interest (i.e., always negative vs. neutral)
P1.posthocBF <-
  data.EEG.trial %>%
  select(participant, cond, P1) %>%
  summarySEwithin(data = ., measurevar = "P1", withinvars = c("participant", "cond")) %>%
  split(.$cond) # split according to condition

for (k in 1:length(scaling.factor)) { # loop through scaling factors

  ### large, dark, negative vs. neutral
  P1.posthocBF.large.dark.negVSneut <- ttestBF(
    P1.posthocBF$`negative large dark`$P1, # first condition
    P1.posthocBF$`neutral large dark`$P1, # second condition
    mu = 0, # null hypothesis (mean difference = 0)
    paired = TRUE, # paired sample test
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.posthocBF.large.dark.negVSneut, # save model
    file = here::here(paste0("EEG/main/models/P1_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### small, dark, negative vs. neutral
  P1.posthocBF.small.dark.negVSneut <- ttestBF(
    P1.posthocBF$`negative small dark`$P1,
    P1.posthocBF$`neutral small dark`$P1,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.posthocBF.small.dark.negVSneut,
    file = here::here(paste0("EEG/main/models/P1_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### large, bright, negative vs. neutral
  P1.posthocBF.large.bright.negVSneut <- ttestBF(
    P1.posthocBF$`negative large bright`$P1,
    P1.posthocBF$`neutral large bright`$P1,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.posthocBF.large.bright.negVSneut,
    file = here::here(paste0("EEG/main/models/P1_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### small, bright, negative vs. neutral
  P1.posthocBF.small.bright.negVSneut <- ttestBF(
    P1.posthocBF$`negative small bright`$P1,
    P1.posthocBF$`neutral small bright`$P1,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.posthocBF.small.bright.negVSneut,
    file = here::here(paste0("EEG/main/models/P1_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds"))
  )
}
# preallocate matrix with all BF10
compare.P1.posthocBF <- matrix(NA, 4, length(scaling.factor))

# preallocate matrix with all % errors
compare.P1.posthocBF.perc.err <- matrix(NA, 4, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {

  ### load models
  P1.posthocBF.large.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  P1.posthocBF.small.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  P1.posthocBF.large.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  P1.posthocBF.small.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### model comparison
  compare.P1.posthocBF[, k] <- c(
    P1.posthocBF.large.dark.negVSneut@bayesFactor$bf,
    P1.posthocBF.small.dark.negVSneut@bayesFactor$bf,
    P1.posthocBF.large.bright.negVSneut@bayesFactor$bf,
    P1.posthocBF.small.bright.negVSneut@bayesFactor$bf
  )

  # percentage of error
  compare.P1.posthocBF.perc.err[, k] <- c(
    P1.posthocBF.large.dark.negVSneut@bayesFactor$error * 100,
    P1.posthocBF.small.dark.negVSneut@bayesFactor$error * 100,
    P1.posthocBF.large.bright.negVSneut@bayesFactor$error * 100,
    P1.posthocBF.small.bright.negVSneut@bayesFactor$error * 100
  )
}

# summary
compare.P1.posthocBF <- data.frame(
  "comparison" = c("large.dark.negVSneut", "small.dark.negVSneut", "large.bright.negVSneut", "small.bright.negVSneut"),
  "nar" = compare.P1.posthocBF[, 1], "nar.p.err" = compare.P1.posthocBF.perc.err[, 1],
  "med" = compare.P1.posthocBF[, 2], "med.p.err" = compare.P1.posthocBF.perc.err[, 2],
  "wid" = compare.P1.posthocBF[, 3], "wid.p.err" = compare.P1.posthocBF.perc.err[, 3]
)

kable(compare.P1.posthocBF)
comparison nar nar.p.err med med.p.err wid wid.p.err
large.dark.negVSneut -1.242476 0.000005 -1.536763 0.000008 -1.852970 0.000009
small.dark.negVSneut -0.083050 0.000004 -0.305802 0.000002 -0.575117 0.000001
large.bright.negVSneut -0.797664 0.000005 -1.063337 0.000005 -1.361115 0.000004
small.bright.negVSneut -1.410897 0.000005 -1.716419 0.000010 -2.039751 0.000011

Follow-up contrasts reveal that the effect of emotion does not explain changes in P1 amplitude:

  • large size, high contrast: 0.215076;
  • small size, high contrast: 0.736532;
  • large size, low contrast: 0.345302;
  • small size, low contrast: 0.179709.

Exploratory model selection

omit.from.full.model nar nar.p.err med med.p.err wid wid.p.err
5 emo 4.48535 39.1146 4.30235 57.3429 4.78873 34.6213
2 cont:emo 2.70675 25.6105 3.54899 54.3118 3.38254 35.0968
1 cont:emo:size 3.32238 26.0123 3.07197 38.4583 3.98027 38.4676
3 emo:size 2.77062 25.2758 2.47795 35.8667 3.17187 31.3004
6 cont -3.45884 37.1962 -3.99157 35.7692 -3.54226 31.1723
4 cont:size -33.87948 27.1131 -34.55479 35.2074 -33.70824 30.3864
7 size -36.90773 28.7470 -36.66896 43.3900 -36.52127 30.2185

Exploratory analyses suggest that omitting the factor emotion from the full model improves fitting by 73.87328 times. Removing the interactions size x contrast x emotion, contrast x emotion, and emotion x size also improves fitting by 34.778268, 21.584342, and 11.91681 times, respectively.
Conversely, omitting the factor contrast or the contrast x size interaction lowers the explanatory value of the resulting model by 54.139968 and 1.01614e+15 times. Finally, omitting the factor size is maximally detrimental, as it would lower the explanatory value of the resulting model by 8.41645e+15 times.

N1

size cont emo N N1 sd se ci
large dark negative 2374 -1.36 3.67 0.08 0.15
large dark neutral 2369 -1.38 3.69 0.08 0.15
large bright negative 2373 -1.03 3.69 0.08 0.15
large bright neutral 2364 -1.13 3.76 0.08 0.15
small dark negative 2383 -1.39 3.67 0.08 0.15
small dark neutral 2376 -1.45 3.71 0.08 0.15
small bright negative 2365 -0.03 3.73 0.08 0.15
small bright neutral 2363 -0.09 3.86 0.08 0.16
for (k in 1:length(scaling.factor)) {

  ### main effects of size and emotion
  N1.sizeplusemo.BF <- lmBF(
    N1 ~ size + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.sizeplusemo.BF,
    file = here::here(paste0("EEG/main/models/N1_sizeplusemo_BF_", scaling.factor[k], ".rds"))
  )

  ### interactive effects of size and emotion
  N1.sizebyemo.BF <- lmBF(
    N1 ~ size * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.sizebyemo.BF,
    file = here::here(paste0("EEG/main/models/N1_sizebyemo_BF_", scaling.factor[k], ".rds"))
  )

  ### main effects of contrast and emotion
  N1.contplusemo.BF <- lmBF(
    N1 ~ cont + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.contplusemo.BF,
    file = here::here(paste0("EEG/main/models/N1_contplusemo_BF_", scaling.factor[k], ".rds"))
  )

  ### interactive effects of contrast and emotion
  N1.contbyemo.BF <- lmBF(
    N1 ~ cont * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.contbyemo.BF,
    file = here::here(paste0("EEG/main/models/N1_contbyemo_BF_", scaling.factor[k], ".rds"))
  )

  ### main effects of size, contrast, and emotion
  N1.sizepluscontplusemo.BF <- lmBF(
    N1 ~ size + cont + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.sizepluscontplusemo.BF,
    file = here::here(paste0("EEG/main/models/N1_sizepluscontplusemo_BF_", scaling.factor[k], ".rds"))
  )

  ### interactive effects of size, contrast, and emotion
  N1.sizebycontbyemo.BF <- lmBF(
    N1 ~ size * cont * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.sizebycontbyemo.BF,
    file = here::here(paste0("EEG/main/models/N1_sizebycontbyemo_BF_", scaling.factor[k], ".rds"))
  )
}
compare.N1.BF <- matrix(NA, 6, length(scaling.factor))
compare.N1.perc.err <- matrix(NA, 6, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) { # loop through scaling factors

  ### load models
  N1.sizeplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_sizeplusemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.sizebyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_sizebyemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.contplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_contplusemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.contbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_contbyemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.sizepluscontplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_sizepluscontplusemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.sizebycontbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_sizebycontbyemo_BF_", scaling.factor[k], ".rds"))
  )

  # BFs
  compare.N1.BF[, k] <- c(
    N1.sizeplusemo.BF@bayesFactor$bf,
    N1.sizebyemo.BF@bayesFactor$bf,
    N1.contplusemo.BF@bayesFactor$bf,
    N1.contbyemo.BF@bayesFactor$bf,
    N1.sizepluscontplusemo.BF@bayesFactor$bf,
    N1.sizebycontbyemo.BF@bayesFactor$bf
  )

  # percentage of error
  compare.N1.perc.err[, k] <- c(
    N1.sizeplusemo.BF@bayesFactor$error * 100,
    N1.sizebyemo.BF@bayesFactor$error * 100,
    N1.contplusemo.BF@bayesFactor$error * 100,
    N1.contbyemo.BF@bayesFactor$error * 100,
    N1.sizepluscontplusemo.BF@bayesFactor$error * 100,
    N1.sizebycontbyemo.BF@bayesFactor$error * 100
  )
}

# summary confirmatory analysis
compare.N1 <- data.frame(
  "model" = c("size + emo", "size x emo", "contr + emo", "cont x emo", "size + cont + emo", "size x cont x emo"),
  "nar" = compare.N1.BF[, 1], "nar.p.err" = compare.N1.perc.err[, 1],
  "med" = compare.N1.BF[, 2], "med.p.err" = compare.N1.perc.err[, 2],
  "wid" = compare.N1.BF[, 3], "wid.p.err" = compare.N1.perc.err[, 3]
)

compare.N1 <- compare.N1[order(compare.N1$med, decreasing = TRUE), ]

kable(compare.N1)
model nar nar.p.err med med.p.err wid wid.p.err
6 size x cont x emo 1136.051 17.90141 1133.523 8.82070 1131.017 8.16727
5 size + cont + emo 1095.365 3.72406 1094.384 4.01831 1093.388 3.95813
3 contr + emo 1054.595 3.24399 1053.978 3.33639 1053.309 3.20027
4 cont x emo 1050.861 3.72420 1049.929 4.09289 1048.905 3.97891
1 size + emo 968.706 3.41983 968.014 3.46315 967.314 3.39620
2 size x emo 964.924 4.17994 963.909 4.63483 962.822 4.64631

When using a JZS prior with scaling factor r = 0.707 placed on standardized effect sizes, the model size x cont x emo ought to be preferred.
The best model (size x cont x emo) explains the observed data 9.95963e+16 times better than the second best model (size + cont + emo).

Paired comparisons

N1.posthocBF <-
  data.EEG.trial %>%
  select(participant, cond, N1) %>%
  summarySEwithin(data = ., measurevar = "N1", withinvars = c("participant", "cond")) %>%
  split(.$cond)

### uncomment to recompute models ###
for (k in 1:length(scaling.factor)) {

  ### large, dark, negative vs. neutral
  N1.posthocBF.large.dark.negVSneut <- ttestBF(
    N1.posthocBF$`negative large dark`$N1,
    N1.posthocBF$`neutral large dark`$N1,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.posthocBF.large.dark.negVSneut, # save model
    file = here::here(paste0("EEG/main/models/N1_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### small, dark, negative vs. neutral
  N1.posthocBF.small.dark.negVSneut <- ttestBF(
    N1.posthocBF$`negative small dark`$N1,
    N1.posthocBF$`neutral small dark`$N1,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.posthocBF.small.dark.negVSneut,
    file = here::here(paste0("EEG/main/models/N1_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### large, bright, negative vs. neutral
  N1.posthocBF.large.bright.negVSneut <- ttestBF(
    N1.posthocBF$`negative large bright`$N1,
    N1.posthocBF$`neutral large bright`$N1,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.posthocBF.large.bright.negVSneut,
    file = here::here(paste0("EEG/main/models/N1_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### small, bright, negative vs. neutral
  N1.posthocBF.small.bright.negVSneut <- ttestBF(
    N1.posthocBF$`negative small bright`$N1,
    N1.posthocBF$`neutral small bright`$N1,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.posthocBF.small.bright.negVSneut,
    file = here::here(paste0("EEG/main/models/N1_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds"))
  )
}
compare.N1.posthocBF <- matrix(NA, 4, length(scaling.factor))
compare.N1.posthocBF.perc.err <- matrix(NA, 4, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {
  N1.posthocBF.large.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  N1.posthocBF.small.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  N1.posthocBF.large.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  N1.posthocBF.small.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  compare.N1.posthocBF[, k] <- c(
    N1.posthocBF.large.dark.negVSneut@bayesFactor$bf,
    N1.posthocBF.small.dark.negVSneut@bayesFactor$bf,
    N1.posthocBF.large.bright.negVSneut@bayesFactor$bf,
    N1.posthocBF.small.bright.negVSneut@bayesFactor$bf
  )

  compare.N1.posthocBF.perc.err[, k] <- c(
    N1.posthocBF.large.dark.negVSneut@bayesFactor$error * 100,
    N1.posthocBF.small.dark.negVSneut@bayesFactor$error * 100,
    N1.posthocBF.large.bright.negVSneut@bayesFactor$error * 100,
    N1.posthocBF.small.bright.negVSneut@bayesFactor$error * 100
  )
}

compare.N1.posthocBF <- data.frame(
  "comparison" = c("large.dark.negVSneut", "small.dark.negVSneut", "large.bright.negVSneut", "small.bright.negVSneut"),
  "nar" = compare.N1.posthocBF[, 1], "nar.p.err" = compare.N1.posthocBF.perc.err[, 1],
  "med" = compare.N1.posthocBF[, 2], "med.p.err" = compare.N1.posthocBF.perc.err[, 2],
  "wid" = compare.N1.posthocBF[, 3], "wid.p.err" = compare.N1.posthocBF.perc.err[, 3]
)

kable(compare.N1.posthocBF)
comparison nar nar.p.err med med.p.err wid wid.p.err
large.dark.negVSneut -1.44404 0.000005 -1.75180 0.000010 -2.07655 0.000012
small.dark.negVSneut -1.34226 0.000005 -1.64318 0.000009 -1.96359 0.000010
large.bright.negVSneut -1.04391 0.000005 -1.32523 0.000006 -1.63315 0.000006
small.bright.negVSneut -1.30534 0.000005 -1.60379 0.000009 -1.92265 0.000010

Follow-up contrasts reveal that the effect of emotion does not explain changes in N1 amplitude:

  • large size, high contrast: 0.173461;
  • small size, high contrast: 0.193365;
  • large size, low contrast: 0.265741;
  • small size, low contrast: 0.201133.

Exploratory model selection

omit.from.full.model nar nar.p.err med med.p.err wid wid.p.err
2 cont:emo 3.73029 11.3004 4.21092 10.08546 4.41023 12.1467
3 emo:size 3.78233 10.9273 4.13584 9.92758 4.57307 22.7891
1 cont:emo:size 3.36395 13.3007 3.99916 9.90268 3.98102 11.7952
5 emo 3.34794 10.9888 3.81375 9.08760 4.15928 13.2988
7 size -41.23282 11.4377 -40.69160 10.02464 -40.54603 12.5904
4 cont:size -51.11622 12.8625 -50.84828 10.26856 -50.72175 12.5625
6 cont -127.35049 11.2607 -126.85935 10.05205 -126.72689 12.1201

Exploratory analyses further suggest that omitting contrast x emotion from the full model improves fitting by 67.418689 times. Removing size x emotion, size x contrast x emotion, and emotion also improves fitting by 62.541803, 54.552148, and 45.320135 times, respectively.
Conversely, omitting the factor size or the contrast x size interaction lowers the explanatory value of the resulting model by 4.70043e+17 and 1.21095e+22 times. Finally, omitting the factor contrast is maximally detrimental, as it would lower the explanatory value of the resulting model by 1.24255e+55 times.

EPN

size cont emo N EPN sd se ci
large dark negative 2374 0.33 3.81 0.08 0.15
large dark neutral 2369 0.53 3.75 0.08 0.15
large bright negative 2373 0.27 3.93 0.08 0.16
large bright neutral 2364 0.58 3.88 0.08 0.16
small dark negative 2383 1.04 3.86 0.08 0.16
small dark neutral 2376 1.37 3.84 0.08 0.15
small bright negative 2365 1.15 4.01 0.08 0.16
small bright neutral 2363 1.09 4.11 0.08 0.17
compare.EPN.BF <- matrix(NA, 6, length(scaling.factor))
compare.EPN.perc.err <- matrix(NA, 6, length(scaling.factor))

for (k in 1:length(scaling.factor)) {

  ### main effects of size and emotion
  EPN.sizeplusemo.BF <- lmBF(
    EPN ~ size + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(EPN.sizeplusemo.BF,
          file = here::here(paste0("EEG/main/models/EPN_sizeplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size and emotion
  EPN.sizebyemo.BF <- lmBF(
    EPN ~ size * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(EPN.sizebyemo.BF,
          file = here::here(paste0("EEG/main/models/EPN_sizebyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of contrast and emotion
  EPN.contplusemo.BF <- lmBF(
    EPN ~ cont + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(EPN.contplusemo.BF,
          file = here::here(paste0("EEG/main/models/EPN_contplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of contrast and emotion
  EPN.contbyemo.BF <- lmBF(
    EPN ~ cont * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(EPN.contbyemo.BF,
          file = here::here(paste0("EEG/main/models/EPN_contbyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of size, contrast, and emotion
  EPN.sizepluscontplusemo.BF <- lmBF(
    EPN ~ size + cont + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(EPN.sizepluscontplusemo.BF,
          file = here::here(paste0("EEG/main/models/EPN_sizepluscontplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size, contrast, and emotion
  EPN.sizebycontbyemo.BF <- lmBF(
    EPN ~ size * cont * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(EPN.sizebycontbyemo.BF,
          file = here::here(paste0("EEG/main/models/EPN_sizebycontbyemo_BF_", scaling.factor[k], ".rds")))

}
compare.EPN.BF <- matrix(NA, 6, length(scaling.factor))
compare.EPN.perc.err <- matrix(NA, 6, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) { # loop through scaling factors

  ### load models
  EPN.sizeplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/EPN_sizeplusemo_BF_", scaling.factor[k], ".rds"))
  )

  EPN.sizebyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/EPN_sizebyemo_BF_", scaling.factor[k], ".rds"))
  )

  EPN.contplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/EPN_contplusemo_BF_", scaling.factor[k], ".rds"))
  )

  EPN.contbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/EPN_contbyemo_BF_", scaling.factor[k], ".rds"))
  )

  EPN.sizepluscontplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/EPN_sizepluscontplusemo_BF_", scaling.factor[k], ".rds"))
  )

  EPN.sizebycontbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/EPN_sizebycontbyemo_BF_", scaling.factor[k], ".rds"))
  )

  # BFs
  compare.EPN.BF[, k] <- c(
    EPN.sizeplusemo.BF@bayesFactor$bf,
    EPN.sizebyemo.BF@bayesFactor$bf,
    EPN.contplusemo.BF@bayesFactor$bf,
    EPN.contbyemo.BF@bayesFactor$bf,
    EPN.sizepluscontplusemo.BF@bayesFactor$bf,
    EPN.sizebycontbyemo.BF@bayesFactor$bf
  )

  # percentage of error
  compare.EPN.perc.err[, k] <- c(
    EPN.sizeplusemo.BF@bayesFactor$error * 100,
    EPN.sizebyemo.BF@bayesFactor$error * 100,
    EPN.contplusemo.BF@bayesFactor$error * 100,
    EPN.contbyemo.BF@bayesFactor$error * 100,
    EPN.sizepluscontplusemo.BF@bayesFactor$error * 100,
    EPN.sizebycontbyemo.BF@bayesFactor$error * 100
  )
  
}

# summary confirmatory analysis
compare.EPN <- data.frame(
  "model" = c("size + emo", "size x emo", "contr + emo", "cont x emo", "size + cont + emo", "size x cont x emo"),
  "nar" = compare.EPN.BF[, 1], "nar.p.err" = compare.EPN.perc.err[, 1],
  "med" = compare.EPN.BF[, 2], "med.p.err" = compare.EPN.perc.err[, 2],
  "wid" = compare.EPN.BF[, 3], "wid.p.err" = compare.EPN.perc.err[, 3]
)

compare.EPN <- compare.EPN[order(compare.EPN$med, decreasing = TRUE), ]

kable(compare.EPN)
model nar nar.p.err med med.p.err wid wid.p.err
1 size + emo 1326.49 2.31222 1325.89 2.40473 1325.18 2.33687
2 size x emo 1323.44 2.87939 1322.39 3.00713 1321.36 2.97128
5 size + cont + emo 1322.75 3.04108 1321.75 3.03328 1320.76 3.06843
6 size x cont x emo 1312.66 6.31408 1310.35 7.02359 1307.86 5.63571
3 contr + emo 1232.33 2.43123 1231.61 2.42752 1230.89 2.44795
4 cont x emo 1229.29 2.93835 1228.24 3.03208 1227.26 3.00893

When using a JZS prior with scaling factor r = 0.707 placed on standardized effect sizes, the model size + emo ought to be preferred.
The best model (size + emo) explains the observed data 33.058994 times better than the second best model (size x emo).

Paired comparisons

EPN.posthocBF <- 
  data.EEG.trial %>%
  select(participant, cond, EPN) %>%
  summarySEwithin(data = ., measurevar = "EPN", withinvars = c("participant", "cond")) %>%
  split(.$cond)

for (k in 1:length(scaling.factor)) {

  ### large, dark, negative vs. neutral
  EPN.posthocBF.large.dark.negVSneut <- ttestBF(
    EPN.posthocBF$`negative large dark`$EPN,
    EPN.posthocBF$`neutral large dark`$EPN,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(EPN.posthocBF.large.dark.negVSneut, # save model
          file = here::here(paste0("EEG/main/models/EPN_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### small, dark, negative vs. neutral
  EPN.posthocBF.small.dark.negVSneut <- ttestBF(
    EPN.posthocBF$`negative small dark`$EPN,
    EPN.posthocBF$`neutral small dark`$EPN,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(EPN.posthocBF.small.dark.negVSneut,
          file = here::here(paste0("EEG/main/models/EPN_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### large, bright, negative vs. neutral
  EPN.posthocBF.large.bright.negVSneut <- ttestBF(
    EPN.posthocBF$`negative large bright`$EPN,
    EPN.posthocBF$`neutral large bright`$EPN,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(EPN.posthocBF.large.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/EPN_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds")))

  ### small, bright, negative vs. neutral
  EPN.posthocBF.small.bright.negVSneut <- ttestBF(
    EPN.posthocBF$`negative small bright`$EPN,
    EPN.posthocBF$`neutral small bright`$EPN,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(EPN.posthocBF.small.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/EPN_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds")))

}
compare.EPN.posthocBF <- matrix(NA, 4, length(scaling.factor))
compare.EPN.posthocBF.perc.err <- matrix(NA, 4, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {
  EPN.posthocBF.large.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/EPN_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  EPN.posthocBF.small.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/EPN_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  EPN.posthocBF.large.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/EPN_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  EPN.posthocBF.small.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/EPN_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  compare.EPN.posthocBF[, k] <- c(
    EPN.posthocBF.large.dark.negVSneut@bayesFactor$bf,
    EPN.posthocBF.small.dark.negVSneut@bayesFactor$bf,
    EPN.posthocBF.large.bright.negVSneut@bayesFactor$bf,
    EPN.posthocBF.small.bright.negVSneut@bayesFactor$bf
  )

  compare.EPN.posthocBF.perc.err[, k] <- c(
    EPN.posthocBF.large.dark.negVSneut@bayesFactor$error * 100,
    EPN.posthocBF.small.dark.negVSneut@bayesFactor$error * 100,
    EPN.posthocBF.large.bright.negVSneut@bayesFactor$error * 100,
    EPN.posthocBF.small.bright.negVSneut@bayesFactor$error * 100
  )
}

compare.EPN.posthocBF <- data.frame(
  "comparison" = c("large.dark.negVSneut", "small.dark.negVSneut", "large.bright.negVSneut", "small.bright.negVSneut"),
  "nar" = compare.EPN.posthocBF[, 1], "nar.p.err" = compare.EPN.posthocBF.perc.err[, 1],
  "med" = compare.EPN.posthocBF[, 2], "med.p.err" = compare.EPN.posthocBF.perc.err[, 2],
  "wid" = compare.EPN.posthocBF[, 3], "wid.p.err" = compare.EPN.posthocBF.perc.err[, 3]
)

kable(compare.EPN.posthocBF)
comparison nar nar.p.err med med.p.err wid wid.p.err
large.dark.negVSneut -0.400217 0.000004 -0.641569 0.000003 -0.923346 0.000002
small.dark.negVSneut 2.410784 0.000000 2.312172 0.000808 2.132030 0.000000
large.bright.negVSneut 1.811765 0.000114 1.686604 0.000000 1.486398 0.000000
small.bright.negVSneut -1.369311 0.000005 -1.672038 0.000009 -1.993603 0.000011

Follow-up contrasts reveal that emotion may explain changes in EPN amplitude when words are presented in small font and high contrast (10.096328), as well as large font and low contrast (5.40111).

Exploratory model selection

omit.from.full.model nar nar.p.err med med.p.err wid wid.p.err
6 cont 3.76437 11.2963 3.95438 10.55308 4.55411 8.35918
4 cont:size 3.48978 11.5879 3.78341 9.74697 4.33640 8.71177
3 emo:size 3.13453 10.5641 3.42454 9.60220 3.90016 8.24739
2 cont:emo 2.91558 11.6178 3.13734 9.36192 3.73162 9.60535
1 cont:emo:size 0.57099 11.4056 0.80486 9.54323 1.18533 7.98465
5 emo -2.91527 10.6151 -2.53273 9.72118 -2.03308 9.34033
7 size -90.39216 11.6568 -90.41695 9.45235 -89.82763 8.50455

Exploratory analyses further suggest that omitting the factor contrast from the full model improves fitting by 52.163543 times. Removing contrast x size, size x emotion, contrast x emotion, and size x contrast x emotion also improves fitting by 43.965611, 30.708517, 23.042374, and 2.236383 times, respectively.
Conversely, omitting the factor emotion lowers the explanatory value of the resulting model by 12.587773 times. Finally, omitting the factor size is maximally detrimental, as it would lower the explanatory value of the resulting model by 1.85175e+39 times.

LPP

size cont emo N LPP sd se ci
large dark negative 2374 0.61 3.04 0.06 0.12
large dark neutral 2369 0.60 2.95 0.06 0.12
large bright negative 2373 0.71 3.03 0.06 0.12
large bright neutral 2364 0.81 3.04 0.06 0.12
small dark negative 2383 0.76 3.09 0.06 0.12
small dark neutral 2376 0.79 3.09 0.06 0.12
small bright negative 2365 1.35 3.31 0.07 0.13
small bright neutral 2363 1.17 3.36 0.07 0.14
for (k in 1:length(scaling.factor)) {

  ### main effects of size and emotion
  LPP.sizeplusemo.BF <- lmBF(
    LPP ~ size + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(LPP.sizeplusemo.BF,
          file = here::here(paste0("EEG/main/models/LPP_sizeplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size and emotion
  LPP.sizebyemo.BF <- lmBF(
    LPP ~ size * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(LPP.sizebyemo.BF,
          file = here::here(paste0("EEG/main/models/LPP_sizebyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of contrast and emotion
  LPP.contplusemo.BF <- lmBF(
    LPP ~ cont + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(LPP.contplusemo.BF,
          file = here::here(paste0("EEG/main/models/LPP_contplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of contrast and emotion
  LPP.contbyemo.BF <- lmBF(
    LPP ~ cont * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(LPP.contbyemo.BF,
          file = here::here(paste0("EEG/main/models/LPP_contbyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of size, contrast, and emotion
  LPP.sizepluscontplusemo.BF <- lmBF(
    LPP ~ size + cont + emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(LPP.sizepluscontplusemo.BF,
          file = here::here(paste0("EEG/main/models/LPP_sizepluscontplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size, contrast, and emotion
  LPP.sizebycontbyemo.BF <- lmBF(
    LPP ~ size * cont * emo + participant,
    data.EEG.trial,
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(LPP.sizebycontbyemo.BF,
          file = here::here(paste0("EEG/main/models/LPP_sizebycontbyemo_BF_", scaling.factor[k], ".rds")))

}
compare.LPP.BF <- matrix(NA, 6, length(scaling.factor))
compare.LPP.perc.err <- matrix(NA, 6, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) { # loop through scaling factors

  ### load models
  LPP.sizeplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/LPP_sizeplusemo_BF_", scaling.factor[k], ".rds"))
  )

  LPP.sizebyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/LPP_sizebyemo_BF_", scaling.factor[k], ".rds"))
  )

  LPP.contplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/LPP_contplusemo_BF_", scaling.factor[k], ".rds"))
  )

  LPP.contbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/LPP_contbyemo_BF_", scaling.factor[k], ".rds"))
  )

  LPP.sizepluscontplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/LPP_sizepluscontplusemo_BF_", scaling.factor[k], ".rds"))
  )

  LPP.sizebycontbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/LPP_sizebycontbyemo_BF_", scaling.factor[k], ".rds"))
  )

  # BFs
  compare.LPP.BF[, k] <- c(
    LPP.sizeplusemo.BF@bayesFactor$bf,
    LPP.sizebyemo.BF@bayesFactor$bf,
    LPP.contplusemo.BF@bayesFactor$bf,
    LPP.contbyemo.BF@bayesFactor$bf,
    LPP.sizepluscontplusemo.BF@bayesFactor$bf,
    LPP.sizebycontbyemo.BF@bayesFactor$bf
  )

  # percentage of error
  compare.LPP.perc.err[, k] <- c(
    LPP.sizeplusemo.BF@bayesFactor$error * 100,
    LPP.sizebyemo.BF@bayesFactor$error * 100,
    LPP.contplusemo.BF@bayesFactor$error * 100,
    LPP.contbyemo.BF@bayesFactor$error * 100,
    LPP.sizepluscontplusemo.BF@bayesFactor$error * 100,
    LPP.sizebycontbyemo.BF@bayesFactor$error * 100
  )
  
}

# summary confirmatory analysis
compare.LPP <- data.frame(
  "model" = c("size + emo", "size x emo", "contr + emo", "cont x emo", "size + cont + emo", "size x cont x emo"),
  "nar" = compare.LPP.BF[, 1], "nar.p.err" = compare.LPP.perc.err[, 1],
  "med" = compare.LPP.BF[, 2], "med.p.err" = compare.LPP.perc.err[, 2],
  "wid" = compare.LPP.BF[, 3], "wid.p.err" = compare.LPP.perc.err[, 3]
)

compare.LPP <- compare.LPP[order(compare.LPP$med, decreasing = TRUE), ]

kable(compare.LPP)
model nar nar.p.err med med.p.err wid wid.p.err
5 size + cont + emo 520.406 11.22651 519.446 17.5976 518.560 11.5647
6 size x cont x emo 516.093 17.08671 513.643 15.5526 511.146 21.9664
1 size + emo 496.161 9.78495 495.544 9.3165 494.769 11.5264
3 contr + emo 494.110 11.09732 493.459 10.7985 492.771 12.4911
2 size x emo 493.421 11.97174 492.545 16.8120 491.177 12.1359
4 cont x emo 490.577 12.05388 489.647 11.9476 488.485 11.5656

When using a JZS prior with scaling factor r = 0.707 placed on standardized effect sizes, the model size + cont + emo ought to be preferred.
The best model (size + cont + emo) explains the observed data 331.122922 times better than the second best model (size x cont x emo).

Paired comparisons

LPP.posthocBF <- 
  data.EEG.trial %>%
  select(participant, cond, LPP) %>%
  summarySEwithin(data = ., measurevar = "LPP", withinvars = c("participant", "cond")) %>%
  split(.$cond)

for (k in 1:length(scaling.factor)) {

  ### large, dark, negative vs. neutral
  LPP.posthocBF.large.dark.negVSneut <- ttestBF(
    LPP.posthocBF$`negative large dark`$LPP,
    LPP.posthocBF$`neutral large dark`$LPP,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(LPP.posthocBF.large.dark.negVSneut, # save model
          file = here::here(paste0("EEG/main/models/LPP_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### small, dark, negative vs. neutral
  LPP.posthocBF.small.dark.negVSneut <- ttestBF(
    LPP.posthocBF$`negative small dark`$LPP,
    LPP.posthocBF$`neutral small dark`$LPP,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(LPP.posthocBF.small.dark.negVSneut,
          file = here::here(paste0("EEG/main/models/LPP_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### large, bright, negative vs. neutral
  LPP.posthocBF.large.bright.negVSneut <- ttestBF(
    LPP.posthocBF$`negative large bright`$LPP,
    LPP.posthocBF$`neutral large bright`$LPP,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(LPP.posthocBF.large.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/LPP_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds")))

  ### small, bright, negative vs. neutral
  LPP.posthocBF.small.bright.negVSneut <- ttestBF(
    LPP.posthocBF$`negative small bright`$LPP,
    LPP.posthocBF$`neutral small bright`$LPP,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(LPP.posthocBF.small.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/LPP_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds")))

}
compare.LPP.posthocBF <- matrix(NA, 4, length(scaling.factor))
compare.LPP.posthocBF.perc.err <- matrix(NA, 4, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {
  LPP.posthocBF.large.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/LPP_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  LPP.posthocBF.small.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/LPP_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  LPP.posthocBF.large.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/LPP_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  LPP.posthocBF.small.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/LPP_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  compare.LPP.posthocBF[, k] <- c(
    LPP.posthocBF.large.dark.negVSneut@bayesFactor$bf,
    LPP.posthocBF.small.dark.negVSneut@bayesFactor$bf,
    LPP.posthocBF.large.bright.negVSneut@bayesFactor$bf,
    LPP.posthocBF.small.bright.negVSneut@bayesFactor$bf
  )

  compare.LPP.posthocBF.perc.err[, k] <- c(
    LPP.posthocBF.large.dark.negVSneut@bayesFactor$error * 100,
    LPP.posthocBF.small.dark.negVSneut@bayesFactor$error * 100,
    LPP.posthocBF.large.bright.negVSneut@bayesFactor$error * 100,
    LPP.posthocBF.small.bright.negVSneut@bayesFactor$error * 100
  )
  
}

compare.LPP.posthocBF <- data.frame(
  "comparison" = c("large.dark.negVSneut", "small.dark.negVSneut", "large.bright.negVSneut", "small.bright.negVSneut"),
  "nar" = compare.LPP.posthocBF[, 1], "nar.p.err" = compare.LPP.posthocBF.perc.err[, 1],
  "med" = compare.LPP.posthocBF[, 2], "med.p.err" = compare.LPP.posthocBF.perc.err[, 2],
  "wid" = compare.LPP.posthocBF[, 3], "wid.p.err" = compare.LPP.posthocBF.perc.err[, 3]
)

kable(compare.LPP.posthocBF)
comparison nar nar.p.err med med.p.err wid wid.p.err
large.dark.negVSneut -1.458985 0.000005 -1.767756 0.000010 -2.093136 0.000012
small.dark.negVSneut -1.388719 0.000005 -1.692748 0.000009 -2.015137 0.000011
large.bright.negVSneut -0.908117 0.000005 -1.180755 0.000005 -1.483057 0.000005
small.bright.negVSneut -0.007331 0.000003 -0.225745 0.000001 -0.492124 0.000001

Follow-up contrasts reveal that the effect of emotion does not explain changes in LPP amplitude:

  • large size, high contrast: 0.170716;
  • small size, high contrast: 0.184013;
  • large size, low contrast: 0.307047;
  • small size, low contrast: 0.797921.

Exploratory model selection

omit.from.full.model nar nar.p.err med med.p.err wid wid.p.err
5 emo 3.97403 28.5906 4.38369 26.1832 5.11209 27.4134
2 cont:emo 3.50601 29.2593 3.75328 28.9674 4.17262 25.7595
3 emo:size 2.69398 30.2241 3.05430 27.2681 3.61779 28.1896
1 cont:emo:size 1.61271 30.2210 1.75300 25.6822 2.70414 30.1085
4 cont:size -3.91812 28.1330 -3.50266 27.6769 -2.90837 26.5374
6 cont -24.71072 25.8800 -24.27263 29.5110 -22.94021 29.9496
7 size -26.34967 31.3473 -25.54122 41.2221 -25.61069 28.5671

Exploratory analyses further suggest that omitting the factor emotion from the full model improves fitting by 80.132853 times. Removing contrast x emotion, size x emotion, and size x contrast x emotion also improves fitting by 42.660828, 21.206335, 5.771892 times, respectively.
Conversely, omitting the contrast x size interaction or the factor contrast lowers the explanatory value of the resulting model by 33.203549 and 3.47912e+10 times, respectively. Finally, omitting the factor size is maximally detrimental, as it would lower the explanatory value of the resulting model by 1.23712e+11 times.



SUPPLEMENTARY ANALYSIS

Visual inspection of the waveforms revealed unexpected latency shifts at the level of P1 and N1 components. This could be a potential source of noise when analyzing mean amplitude values. Therefore, we ran the same analyses as above but using, as dependent variables, peak amplitude and latency values of P1 and N1.

# load data
data.P1.peakamp <- read.table(here::here("EEG/main/data/P1_locpeakamp.txt"), header = TRUE)
data.P1.peaklat <- read.table(here::here("EEG/main/data/P1_locpeaklat.txt"), header = TRUE)
data.N1.peakamp <- read.table(here::here("EEG/main/data/N1_locpeakamp.txt"), header = TRUE)
data.N1.peaklat <- read.table(here::here("EEG/main/data/N1_locpeaklat.txt"), header = TRUE)

data.peaks <-
  rbind(data.P1.peakamp, data.P1.peaklat, data.N1.peakamp, data.N1.peaklat) %>%
  mutate(
    var = rep(c("peak.amp", "peak.lat", "peak.amp", "peak.lat"), each = 320), # peak amplitude or latency?
    bini = factor(bini), # convert condition numbers as factors
    size = revalue(
      bini, # main effect of size
      c(
        "1" = "large", "2" = "small", "3" = "large", "4" = "small",
        "5" = "large", "6" = "small", "7" = "large", "8" = "small"
      )
    ),
    size = relevel(size, ref = "large"), # re-reference size to large
    cont = revalue(
      bini, # main effect of contrast
      c(
        "1" = "dark", "2" = "dark", "3" = "bright", "4" = "bright",
        "5" = "dark", "6" = "dark", "7" = "bright", "8" = "bright"
      )
    ),
    cont = relevel(cont, ref = "dark"), # re-reference contrast to dark
    emo = revalue(
      bini, # main effect of emotion
      c(
        "1" = "negative", "2" = "negative", "3" = "negative", "4" = "negative",
        "5" = "neutral", "6" = "neutral", "7" = "neutral", "8" = "neutral"
      )
    ),
    emo = relevel(emo, ref = "neutral") # re-reference emotion to neutral
  ) %>%
  unite(visualfeats, c("size", "cont"), sep = " ", remove = FALSE) %>% # for graph
  mutate(visualfeats = as.factor(ifelse(visualfeats == "NA NA", "localizer", visualfeats))) %>%
  # select, reorder, and rename variables
  dplyr::select(participant = ERPset, component = chlabel, var, condition = binlabel, size, cont, emo, visualfeats, value)

rm(data.P1.peakamp, data.P1.peaklat, data.N1.peakamp, data.N1.peaklat) # remove unnecessary data frames

P1

Peak amplitude

Peak amplitude: positive value larger than the 3 points (~10 ms at 256 Hz sampling rate) on either side of the peak.
Time window: 66-148 ms post-stimulus onset.
Electrode cluster: P7 P9 PO7 O1 O2 PO8 P8 P10.

size cont emo visualfeats N value sd se ci
large dark neutral large dark 40 3.35137 1.025301 0.162114 0.327907
large dark negative large dark 40 3.29277 1.064229 0.168269 0.340357
large bright neutral large bright 40 3.00358 0.940474 0.148702 0.300778
large bright negative large bright 40 3.13585 0.873157 0.138058 0.279249
small dark neutral small dark 40 2.66942 0.941481 0.148861 0.301100
small dark negative small dark 40 2.90865 0.833393 0.131771 0.266532
small bright neutral small bright 40 1.29638 1.664701 0.263212 0.532397
small bright negative small bright 40 1.46732 1.222649 0.193318 0.391022

for (k in 1:length(scaling.factor)) {

  ### main effects of size and emotion
  P1.peakamp.sizeplusemo.BF <- lmBF(
    value ~ size + emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peakamp.sizeplusemo.BF, # save model
          file = here::here(paste0("EEG/main/models/P1_peakamp_sizeplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size and emotion
  P1.peakamp.sizebyemo.BF <- lmBF(
    value ~ size * emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peakamp.sizebyemo.BF,
          file = here::here(paste0("EEG/main/models/P1_peakamp_sizebyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of contrast and emotion
  P1.peakamp.contplusemo.BF <- lmBF(
    value ~ cont + emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peakamp.contplusemo.BF,
          file = here::here(paste0("EEG/main/models/P1_peakamp_contplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of contrast and emotion
  P1.peakamp.contbyemo.BF <- lmBF(
    value ~ cont * emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peakamp.contbyemo.BF,
          file = here::here(paste0("EEG/main/models/P1_peakamp_contbyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of size, contrast, and emotion
  P1.peakamp.sizepluscontplusemo.BF <- lmBF(
    value ~ size + cont + emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peakamp.sizepluscontplusemo.BF,
          file = here::here(paste0("EEG/main/models/P1_peakamp_sizepluscontplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size, contrast, and emotion
  P1.peakamp.sizebycontbyemo.BF <- lmBF(
    value ~ size * cont * emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peakamp.sizebycontbyemo.BF,
          file = here::here(paste0("EEG/main/models/P1_peakamp_sizebycontbyemo_BF_", scaling.factor[k], ".rds")))

}
compare.P1.peakamp.BF <- matrix(NA, 6, length(scaling.factor))
compare.P1.peakamp.perc.err <- matrix(NA, 6, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {

  ### load models
  P1.peakamp.sizeplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peakamp_sizeplusemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.peakamp.sizebyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peakamp_sizebyemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.peakamp.contplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peakamp_contplusemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.peakamp.contbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peakamp_contbyemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.peakamp.sizepluscontplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peakamp_sizepluscontplusemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.peakamp.sizebycontbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peakamp_sizebycontbyemo_BF_", scaling.factor[k], ".rds"))
  )

  compare.P1.peakamp.BF[, k] <- c(
    P1.peakamp.sizeplusemo.BF@bayesFactor$bf,
    P1.peakamp.sizebyemo.BF@bayesFactor$bf,
    P1.peakamp.contplusemo.BF@bayesFactor$bf,
    P1.peakamp.contbyemo.BF@bayesFactor$bf,
    P1.peakamp.sizepluscontplusemo.BF@bayesFactor$bf,
    P1.peakamp.sizebycontbyemo.BF@bayesFactor$bf
  )

  compare.P1.peakamp.perc.err[, k] <- c(
    P1.peakamp.sizeplusemo.BF@bayesFactor$error * 100,
    P1.peakamp.sizebyemo.BF@bayesFactor$error * 100,
    P1.peakamp.contplusemo.BF@bayesFactor$error * 100,
    P1.peakamp.contbyemo.BF@bayesFactor$error * 100,
    P1.peakamp.sizepluscontplusemo.BF@bayesFactor$error * 100,
    P1.peakamp.sizebycontbyemo.BF@bayesFactor$error * 100
  )
  
}

compare.P1.peakamp <- data.frame(
  "model" = c("size + emo", "size x emo", "contr + emo", "cont x emo", "size + cont + emo", "size x cont x emo"),
  "nar" = compare.P1.peakamp.BF[, 1], "nar.p.err" = compare.P1.peakamp.perc.err[, 1],
  "med" = compare.P1.peakamp.BF[, 2], "med.p.err" = compare.P1.peakamp.perc.err[, 2],
  "wid" = compare.P1.peakamp.BF[, 3], "wid.p.err" = compare.P1.peakamp.perc.err[, 3]
)

compare.P1.peakamp <- compare.P1.peakamp[order(compare.P1.peakamp$med, decreasing = TRUE), ]

kable(compare.P1.peakamp)
model nar nar.p.err med med.p.err wid wid.p.err
6 size x cont x emo 113.0472 2.082444 111.6635 2.226273 109.9492 2.505049
5 size + cont + emo 110.3954 0.972598 109.9963 0.996991 109.4234 1.068096
1 size + emo 94.8978 0.864897 94.5627 0.880644 94.1140 0.904550
2 size x emo 93.2742 1.068691 92.6474 1.096823 91.8743 1.137737
3 contr + emo 82.2755 0.853039 81.8464 0.885007 81.3105 0.927919
4 cont x emo 80.5544 1.060617 79.7795 1.120119 78.9037 1.169013

When using a JZS prior with scaling factor r = 0.707 placed on standardized effect sizes, the model size x cont x emo ought to be preferred.
The best model (size x cont x emo) explains the observed data 5.29724 times better than the second best model (size + cont + emo).

Paired comparisons

P1.peakamp.posthocBF <- 
  subset(data.peaks, component == "P1" & var == "peak.amp") %>%
  select(participant, condition, value) %>%
  split(.$condition)

for (k in 1:length(scaling.factor)) {

  ### large, dark, negative vs. neutral
  P1.peakamp.posthocBF.large.dark.negVSneut <- ttestBF(
    P1.peakamp.posthocBF$`negativeLargeDark`$value,
    P1.peakamp.posthocBF$`neutralLargeDark`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.peakamp.posthocBF.large.dark.negVSneut, # save model
          file = here::here(paste0("EEG/main/models/P1_peakamp_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### small, dark, negative vs. neutral
  P1.peakamp.posthocBF.small.dark.negVSneut <- ttestBF(
    P1.peakamp.posthocBF$`negativeSmallDark`$value,
    P1.peakamp.posthocBF$`neutralSmallDark`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.peakamp.posthocBF.small.dark.negVSneut,
          file = here::here(paste0("EEG/main/models/P1_peakamp_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### large, bright, negative vs. neutral
  P1.peakamp.posthocBF.large.bright.negVSneut <- ttestBF(
    P1.peakamp.posthocBF$`negativeLargeBright`$value,
    P1.peakamp.posthocBF$`neutralLargeBright`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.peakamp.posthocBF.large.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/P1_peakamp_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds")))

  ### small, bright, negative vs. neutral
  P1.peakamp.posthocBF.small.bright.negVSneut <- ttestBF(
    P1.peakamp.posthocBF$`negativeSmallBright`$value,
    P1.peakamp.posthocBF$`neutralSmallBright`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.peakamp.posthocBF.small.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/P1_peakamp_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds")))

}
compare.P1.peakamp.posthocBF <- matrix(NA, 4, length(scaling.factor))
compare.P1.peakamp.posthocBF.perc.err <- matrix(NA, 4, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {

  ### load models
  P1.peakamp.posthocBF.large.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_peakamp_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  P1.peakamp.posthocBF.small.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_peakamp_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  P1.peakamp.posthocBF.large.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_peakamp_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  P1.peakamp.posthocBF.small.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_peakamp_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### model comparison
  compare.P1.peakamp.posthocBF[, k] <- c(
    P1.peakamp.posthocBF.large.dark.negVSneut@bayesFactor$bf,
    P1.peakamp.posthocBF.small.dark.negVSneut@bayesFactor$bf,
    P1.peakamp.posthocBF.large.bright.negVSneut@bayesFactor$bf,
    P1.peakamp.posthocBF.small.bright.negVSneut@bayesFactor$bf
  )

  # percentage of error
  compare.P1.peakamp.posthocBF.perc.err[, k] <- c(
    P1.peakamp.posthocBF.large.dark.negVSneut@bayesFactor$error * 100,
    P1.peakamp.posthocBF.small.dark.negVSneut@bayesFactor$error * 100,
    P1.peakamp.posthocBF.large.bright.negVSneut@bayesFactor$error * 100,
    P1.peakamp.posthocBF.small.bright.negVSneut@bayesFactor$error * 100
  )
}

# summary
compare.P1.peakamp.posthocBF <- data.frame(
  "comparison" = c("large.dark.negVSneut", "small.dark.negVSneut", "large.bright.negVSneut", "small.bright.negVSneut"),
  "nar" = compare.P1.peakamp.posthocBF[, 1], "nar.p.err" = compare.P1.peakamp.posthocBF.perc.err[, 1],
  "med" = compare.P1.peakamp.posthocBF[, 2], "med.p.err" = compare.P1.peakamp.posthocBF.perc.err[, 2],
  "wid" = compare.P1.peakamp.posthocBF[, 3], "wid.p.err" = compare.P1.peakamp.posthocBF.perc.err[, 3]
)

kable(compare.P1.peakamp.posthocBF)
comparison nar nar.p.err med med.p.err wid wid.p.err
large.dark.negVSneut -1.32816 0.000005 -1.628132 0.000009 -1.947953 0.000010
small.dark.negVSneut -0.07760 0.000004 -0.300039 0.000002 -0.569142 0.000001
large.bright.negVSneut -1.22201 0.000005 -1.514945 0.000008 -1.830291 0.000008
small.bright.negVSneut -1.19346 0.000005 -1.484519 0.000008 -1.798667 0.000008

Follow-up contrasts reveal that the effect of emotion does not explain changes in P1 peak amplitude:

  • large size, high contrast: 0.196296;
  • small size, high contrast: 0.74079;
  • large size, low contrast: 0.21982;
  • small size, low contrast: 0.226611.

Exploratory model selection

omit.from.full.model nar nar.p.err med med.p.err wid wid.p.err
2 cont:emo 1.72525 2.53469 2.01233 2.84692 2.39872 3.21703
5 emo 1.65298 2.53784 1.92999 2.88048 2.34748 3.24195
3 emo:size 1.58722 2.52966 1.81079 2.83664 2.20254 3.25555
1 cont:emo:size 1.35922 2.54836 1.60462 2.87772 1.96946 3.29696
4 cont:size -7.33879 2.60862 -7.36829 2.93751 -7.16748 3.24773
6 cont -16.67293 2.65234 -16.72574 2.93242 -16.59378 3.23379
7 size -30.03498 2.65536 -30.21739 3.00605 -30.15072 3.38574

Exploratory analyses suggest that omitting contrast x emotion from the full model improves fitting by 7.480755 times. Removing emotion, size x emotion, and size x contrast x emotion also improves fitting by 6.889455, 6.115285, and 4.975987 times, respectively.
Conversely, omitting contrast x size or contrast lowers the explanatory value of the resulting model by 1584.929377 and 1.8361e+07 times, respectively. Finally, omitting the factor size is maximally detrimental, as it would lower the explanatory value of the resulting model by 1.32815e+13 times.

Peak latency

Peak latency: time (in ms) of the positive value larger than the 3 points (~10 ms at 256 Hz sampling rate) on either side of the peak.
Time window: 66-148 ms post-stimulus onset.
Electrode cluster: P7 P9 PO7 O1 O2 PO8 P8 P10.

size cont emo visualfeats N value sd se ci
large dark neutral large dark 40 103.711 13.59364 2.14934 4.34746
large dark negative large dark 40 102.344 15.28566 2.41687 4.88859
large bright neutral large bright 40 123.633 21.13074 3.34106 6.75794
large bright negative large bright 40 126.270 17.38954 2.74953 5.56145
small dark neutral small dark 40 113.184 9.92574 1.56940 3.17441
small dark negative small dark 40 116.699 13.39870 2.11852 4.28511
small bright neutral small bright 40 115.430 27.00762 4.27028 8.63745
small bright negative small bright 40 116.992 22.47403 3.55346 7.18754

for (k in 1:length(scaling.factor)) {

  ### main effects of size and emotion
  P1.peaklat.sizeplusemo.BF <- lmBF(
    value ~ size + emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peaklat.sizeplusemo.BF, # save model
          file = here::here(paste0("EEG/main/models/P1_peaklat_sizeplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size and emotion
  P1.peaklat.sizebyemo.BF <- lmBF(
    value ~ size * emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peaklat.sizebyemo.BF,
          file = here::here(paste0("EEG/main/models/P1_peaklat_sizebyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of contrast and emotion
  P1.peaklat.contplusemo.BF <- lmBF(
    value ~ cont + emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peaklat.contplusemo.BF,
          file = here::here(paste0("EEG/main/models/P1_peaklat_contplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of contrast and emotion
  P1.peaklat.contbyemo.BF <- lmBF(
    value ~ cont * emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peaklat.contbyemo.BF,
          file = here::here(paste0("EEG/main/models/P1_peaklat_contbyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of size, contrast, and emotion
  P1.peaklat.sizepluscontplusemo.BF <- lmBF(
    value ~ size + cont + emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peaklat.sizepluscontplusemo.BF,
          file = here::here(paste0("EEG/main/models/P1_peaklat_sizepluscontplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size, contrast, and emotion
  P1.peaklat.sizebycontbyemo.BF <- lmBF(
    value ~ size * cont * emo + participant,
    subset(data.peaks, component == "P1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(P1.peaklat.sizebycontbyemo.BF,
          file = here::here(paste0("EEG/main/models/P1_peaklat_sizebycontbyemo_BF_", scaling.factor[k], ".rds")))

}
compare.P1.peaklat.BF <- matrix(NA, 6, length(scaling.factor))
compare.P1.peaklat.perc.err <- matrix(NA, 6, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {

  ### load models
  P1.peaklat.sizeplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peaklat_sizeplusemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.peaklat.sizebyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peaklat_sizebyemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.peaklat.contplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peaklat_contplusemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.peaklat.contbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peaklat_contbyemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.peaklat.sizepluscontplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peaklat_sizepluscontplusemo_BF_", scaling.factor[k], ".rds"))
  )

  P1.peaklat.sizebycontbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/P1_peaklat_sizebycontbyemo_BF_", scaling.factor[k], ".rds"))
  )

  compare.P1.peaklat.BF[, k] <- c(
    P1.peaklat.sizeplusemo.BF@bayesFactor$bf,
    P1.peaklat.sizebyemo.BF@bayesFactor$bf,
    P1.peaklat.contplusemo.BF@bayesFactor$bf,
    P1.peaklat.contbyemo.BF@bayesFactor$bf,
    P1.peaklat.sizepluscontplusemo.BF@bayesFactor$bf,
    P1.peaklat.sizebycontbyemo.BF@bayesFactor$bf
  )

  compare.P1.peaklat.perc.err[, k] <- c(
    P1.peaklat.sizeplusemo.BF@bayesFactor$error * 100,
    P1.peaklat.sizebyemo.BF@bayesFactor$error * 100,
    P1.peaklat.contplusemo.BF@bayesFactor$error * 100,
    P1.peaklat.contbyemo.BF@bayesFactor$error * 100,
    P1.peaklat.sizepluscontplusemo.BF@bayesFactor$error * 100,
    P1.peaklat.sizebycontbyemo.BF@bayesFactor$error * 100
  )
}

compare.P1.peaklat <- data.frame(
  "model" = c("size + emo", "size x emo", "contr + emo", "cont x emo", "size + cont + emo", "size x cont x emo"),
  "nar" = compare.P1.peaklat.BF[, 1], "nar.p.err" = compare.P1.peaklat.perc.err[, 1],
  "med" = compare.P1.peaklat.BF[, 2], "med.p.err" = compare.P1.peaklat.perc.err[, 2],
  "wid" = compare.P1.peaklat.BF[, 3], "wid.p.err" = compare.P1.peaklat.perc.err[, 3]
)

compare.P1.peaklat <- compare.P1.peaklat[order(compare.P1.peaklat$med, decreasing = TRUE), ]

kable(compare.P1.peaklat)
model nar nar.p.err med med.p.err wid wid.p.err
6 size x cont x emo 10.39443 3.93015 8.64674 4.27277 6.61916 4.98297
3 contr + emo 7.84221 1.90041 7.38171 1.95927 6.84259 2.03116
4 cont x emo 6.09696 2.25124 5.32543 2.44768 4.40629 2.55865
5 size + cont + emo 6.00173 2.28447 5.21400 2.38414 4.36215 2.59607
1 size + emo -4.87032 2.48111 -5.49833 2.45132 -6.18859 2.48785
2 size x emo -6.61149 2.99039 -7.54807 3.01909 -8.55482 3.00385

When using a JZS prior with scaling factor r = 0.707 placed on standardized effect sizes, the model size x cont x emo ought to be preferred.
The best model (size x cont x emo) explains the observed data 3.543206 times better than the second best model (contr + emo).

Paired comparisons

P1.peaklat.posthocBF <- 
  subset(data.peaks, component == "P1" & var == "peak.lat") %>%
  select(participant, condition, value) %>%
  split(.$condition)

for (k in 1:length(scaling.factor)) {

  ### large, dark, negative vs. neutral
  P1.peaklat.posthocBF.large.dark.negVSneut <- ttestBF(
    P1.peaklat.posthocBF$`negativeLargeDark`$value,
    P1.peaklat.posthocBF$`neutralLargeDark`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.peaklat.posthocBF.large.dark.negVSneut, # save model
          file = here::here(paste0("EEG/main/models/P1_peaklat_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### small, dark, negative vs. neutral
  P1.peaklat.posthocBF.small.dark.negVSneut <- ttestBF(
    P1.peaklat.posthocBF$`negativeSmallDark`$value,
    P1.peaklat.posthocBF$`neutralSmallDark`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.peaklat.posthocBF.small.dark.negVSneut,
          file = here::here(paste0("EEG/main/models/P1_peaklat_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### large, bright, negative vs. neutral
  P1.peaklat.posthocBF.large.bright.negVSneut <- ttestBF(
    P1.peaklat.posthocBF$`negativeLargeBright`$value,
    P1.peaklat.posthocBF$`neutralLargeBright`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.peaklat.posthocBF.large.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/P1_peaklat_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds")))

  ### small, bright, negative vs. neutral
  P1.peaklat.posthocBF.small.bright.negVSneut <- ttestBF(
    P1.peaklat.posthocBF$`negativeSmallBright`$value,
    P1.peaklat.posthocBF$`neutralSmallBright`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(P1.peaklat.posthocBF.small.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/P1_peaklat_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds")))

}
compare.P1.peaklat.posthocBF <- matrix(NA, 4, length(scaling.factor))
compare.P1.peaklat.posthocBF.perc.err <- matrix(NA, 4, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {

  ### load models
  P1.peaklat.posthocBF.large.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_peaklat_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  P1.peaklat.posthocBF.small.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_peaklat_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  P1.peaklat.posthocBF.large.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_peaklat_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  P1.peaklat.posthocBF.small.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/P1_peaklat_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### model comparison
  compare.P1.peaklat.posthocBF[, k] <- c(
    P1.peaklat.posthocBF.large.dark.negVSneut@bayesFactor$bf,
    P1.peaklat.posthocBF.small.dark.negVSneut@bayesFactor$bf,
    P1.peaklat.posthocBF.large.bright.negVSneut@bayesFactor$bf,
    P1.peaklat.posthocBF.small.bright.negVSneut@bayesFactor$bf
  )

  # percentage of error
  compare.P1.peaklat.posthocBF.perc.err[, k] <- c(
    P1.peaklat.posthocBF.large.dark.negVSneut@bayesFactor$error * 100,
    P1.peaklat.posthocBF.small.dark.negVSneut@bayesFactor$error * 100,
    P1.peaklat.posthocBF.large.bright.negVSneut@bayesFactor$error * 100,
    P1.peaklat.posthocBF.small.bright.negVSneut@bayesFactor$error * 100
  )
}

# summary
compare.P1.peaklat.posthocBF <- data.frame(
  "comparison" = c("large.dark.negVSneut", "small.dark.negVSneut", "large.bright.negVSneut", "small.bright.negVSneut"),
  "nar" = compare.P1.peaklat.posthocBF[, 1], "nar.p.err" = compare.P1.peaklat.posthocBF.perc.err[, 1],
  "med" = compare.P1.peaklat.posthocBF[, 2], "med.p.err" = compare.P1.peaklat.posthocBF.perc.err[, 2],
  "wid" = compare.P1.peaklat.posthocBF[, 3], "wid.p.err" = compare.P1.peaklat.posthocBF.perc.err[, 3]
)

kable(compare.P1.peaklat.posthocBF)
comparison nar nar.p.err med med.p.err wid wid.p.err
large.dark.negVSneut -1.324323 0.000005 -1.624043 0.000009 -1.943702 0.000010
small.dark.negVSneut -0.436943 0.000004 -0.680494 0.000003 -0.963732 0.000002
large.bright.negVSneut -1.288675 0.000005 -1.586022 0.000008 -1.904176 0.000009
small.bright.negVSneut -1.432421 0.000005 -1.739395 0.000010 -2.063643 0.000012

Follow-up contrasts reveal that the effect of emotion does not explain changes in P1 peak latency:

  • large size, high contrast: 0.1971;
  • small size, high contrast: 0.506367;
  • large size, low contrast: 0.204738;
  • small size, low contrast: 0.175627.

Exploratory model selection

omit.from.full.model nar nar.p.err med med.p.err wid wid.p.err
5 emo 1.83285 5.23975 2.15809 6.90986 2.43201 6.27148
7 size 1.82995 5.12262 2.12552 6.99069 2.47590 6.37073
3 emo:size 1.63757 4.98234 2.00774 6.91695 2.32132 6.39586
2 cont:emo 1.71518 5.08172 1.97646 6.72197 2.33141 6.17120
1 cont:emo:size 1.24054 5.14424 1.45722 6.87122 1.75279 6.16118
4 cont:size -9.05984 5.82549 -9.19382 6.99118 -8.97445 6.56233
6 cont -11.88366 5.75009 -11.77912 7.21225 -11.63297 6.35138

Exploratory analyses suggest that omitting emotion from the full model improves fitting by 8.654611 times. Removing size, contrast x emotion, size x emotion, and size x contrast x emotion also improves fitting by 8.377287, 7.446467, 7.217118, and 4.294014 times, respectively.
Conversely, omitting contrast x size lowers the explanatory value of the resulting model by 9836.163801 times. Finally, omitting the factor contrast is maximally detrimental, as it would lower the explanatory value of the resulting model by 1.30499e+05 times.

N1

Peak amplitude

Peak amplitude: negative value larger than the 3 points (~10 ms at 256 Hz sampling rate) on either side of the peak.
Time window: 150-260 ms post-stimulus onset.
Electrode cluster: TP7 P7 P9 TP8 P8 P10.

size cont emo visualfeats N value sd se ci
large dark neutral large dark 40 -3.497500 1.130311 0.178718 0.361491
large dark negative large dark 40 -3.506750 1.046491 0.165465 0.334684
large bright neutral large bright 40 -3.051325 0.835515 0.132107 0.267211
large bright negative large bright 40 -3.067950 0.711714 0.112532 0.227617
small dark neutral small dark 40 -3.302025 1.089445 0.172256 0.348421
small dark negative small dark 40 -3.243100 0.881943 0.139447 0.282059
small bright neutral small bright 40 -1.244675 1.309893 0.207112 0.418924
small bright negative small bright 40 -0.992025 1.719284 0.271843 0.549854

for (k in 1:length(scaling.factor)) {

  ### main effects of size and emotion
  N1.peakamp.sizeplusemo.BF <- lmBF(
    value ~ size + emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peakamp.sizeplusemo.BF, # save model
          file = here::here(paste0("EEG/main/models/N1_peakamp_sizeplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size and emotion
  N1.peakamp.sizebyemo.BF <- lmBF(
    value ~ size * emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peakamp.sizebyemo.BF,
          file = here::here(paste0("EEG/main/models/N1_peakamp_sizebyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of contrast and emotion
  N1.peakamp.contplusemo.BF <- lmBF(
    value ~ cont + emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peakamp.contplusemo.BF,
          file = here::here(paste0("EEG/main/models/N1_peakamp_contplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of contrast and emotion
  N1.peakamp.contbyemo.BF <- lmBF(
    value ~ cont * emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peakamp.contbyemo.BF,
          file = here::here(paste0("EEG/main/models/N1_peakamp_contbyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of size, contrast, and emotion
  N1.peakamp.sizepluscontplusemo.BF <- lmBF(
    value ~ size + cont + emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peakamp.sizepluscontplusemo.BF,
          file = here::here(paste0("EEG/main/models/N1_peakamp_sizepluscontplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size, contrast, and emotion
  N1.peakamp.sizebycontbyemo.BF <- lmBF(
    value ~ size * cont * emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.amp"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peakamp.sizebycontbyemo.BF,
          file = here::here(paste0("EEG/main/models/N1_peakamp_sizebycontbyemo_BF_", scaling.factor[k], ".rds")))

}
compare.N1.peakamp.BF <- matrix(NA, 6, length(scaling.factor))
compare.N1.peakamp.perc.err <- matrix(NA, 6, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {

  ### load models
  N1.peakamp.sizeplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peakamp_sizeplusemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.peakamp.sizebyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peakamp_sizebyemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.peakamp.contplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peakamp_contplusemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.peakamp.contbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peakamp_contbyemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.peakamp.sizepluscontplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peakamp_sizepluscontplusemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.peakamp.sizebycontbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peakamp_sizebycontbyemo_BF_", scaling.factor[k], ".rds"))
  )

  compare.N1.peakamp.BF[, k] <- c(
    N1.peakamp.sizeplusemo.BF@bayesFactor$bf,
    N1.peakamp.sizebyemo.BF@bayesFactor$bf,
    N1.peakamp.contplusemo.BF@bayesFactor$bf,
    N1.peakamp.contbyemo.BF@bayesFactor$bf,
    N1.peakamp.sizepluscontplusemo.BF@bayesFactor$bf,
    N1.peakamp.sizebycontbyemo.BF@bayesFactor$bf
  )

  compare.N1.peakamp.perc.err[, k] <- c(
    N1.peakamp.sizeplusemo.BF@bayesFactor$error * 100,
    N1.peakamp.sizebyemo.BF@bayesFactor$error * 100,
    N1.peakamp.contplusemo.BF@bayesFactor$error * 100,
    N1.peakamp.contbyemo.BF@bayesFactor$error * 100,
    N1.peakamp.sizepluscontplusemo.BF@bayesFactor$error * 100,
    N1.peakamp.sizebycontbyemo.BF@bayesFactor$error * 100
  )
}

compare.N1.peakamp <- data.frame(
  "model" = c("size + emo", "size x emo", "contr + emo", "cont x emo", "size + cont + emo", "size x cont x emo"),
  "nar" = compare.N1.peakamp.BF[, 1], "nar.p.err" = compare.N1.peakamp.perc.err[, 1],
  "med" = compare.N1.peakamp.BF[, 2], "med.p.err" = compare.N1.peakamp.perc.err[, 2],
  "wid" = compare.N1.peakamp.BF[, 3], "wid.p.err" = compare.N1.peakamp.perc.err[, 3]
)

compare.N1.peakamp <- compare.N1.peakamp[order(compare.N1.peakamp$med, decreasing = TRUE), ]

kable(compare.N1.peakamp)
model nar nar.p.err med med.p.err wid wid.p.err
6 size x cont x emo 132.5116 2.188652 131.3568 2.215321 129.8644 2.308238
5 size + cont + emo 120.4792 1.014891 120.2075 1.004733 119.7191 1.040088
3 contr + emo 96.7546 0.868479 96.4366 0.873540 96.0086 0.889878
4 cont x emo 95.0109 1.083816 94.4108 1.126364 93.6376 1.141639
1 size + emo 86.8748 0.852969 86.4728 0.874806 85.9936 0.916449
2 size x emo 85.2270 1.087511 84.5178 1.114712 83.7137 1.168779

When using a JZS prior with scaling factor r = 0.707 placed on standardized effect sizes, the model size x cont x emo ought to be preferred.
The best model (size x cont x emo) explains the observed data 6.95162e+04 times better than the second best model (size + cont + emo).

Paired comparisons

N1.peakamp.posthocBF <- 
  subset(data.peaks, component == "N1" & var == "peak.amp") %>%
  select(participant, condition, value) %>%
  split(.$condition)

for (k in 1:length(scaling.factor)) {

  ### large, dark, negative vs. neutral
  N1.peakamp.posthocBF.large.dark.negVSneut <- ttestBF(
    N1.peakamp.posthocBF$`negativeLargeDark`$value,
    N1.peakamp.posthocBF$`neutralLargeDark`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.peakamp.posthocBF.large.dark.negVSneut, # save model
          file = here::here(paste0("EEG/main/models/N1_peakamp_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### small, dark, negative vs. neutral
  N1.peakamp.posthocBF.small.dark.negVSneut <- ttestBF(
    N1.peakamp.posthocBF$`negativeSmallDark`$value,
    N1.peakamp.posthocBF$`neutralSmallDark`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.peakamp.posthocBF.small.dark.negVSneut,
          file = here::here(paste0("EEG/main/models/N1_peakamp_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### large, bright, negative vs. neutral
  N1.peakamp.posthocBF.large.bright.negVSneut <- ttestBF(
    N1.peakamp.posthocBF$`negativeLargeBright`$value,
    N1.peakamp.posthocBF$`neutralLargeBright`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.peakamp.posthocBF.large.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/N1_peakamp_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds")))

  ### small, bright, negative vs. neutral
  N1.peakamp.posthocBF.small.bright.negVSneut <- ttestBF(
    N1.peakamp.posthocBF$`negativeSmallBright`$value,
    N1.peakamp.posthocBF$`neutralSmallBright`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.peakamp.posthocBF.small.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/N1_peakamp_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds")))

}
compare.N1.peakamp.posthocBF <- matrix(NA, 4, length(scaling.factor))
compare.N1.peakamp.posthocBF.perc.err <- matrix(NA, 4, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {

  ### load models
  N1.peakamp.posthocBF.large.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_peakamp_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  N1.peakamp.posthocBF.small.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_peakamp_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  N1.peakamp.posthocBF.large.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_peakamp_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  N1.peakamp.posthocBF.small.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_peakamp_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### model comparison
  compare.N1.peakamp.posthocBF[, k] <- c(
    N1.peakamp.posthocBF.large.dark.negVSneut@bayesFactor$bf,
    N1.peakamp.posthocBF.small.dark.negVSneut@bayesFactor$bf,
    N1.peakamp.posthocBF.large.bright.negVSneut@bayesFactor$bf,
    N1.peakamp.posthocBF.small.bright.negVSneut@bayesFactor$bf
  )

  # percentage of error
  compare.N1.peakamp.posthocBF.perc.err[, k] <- c(
    N1.peakamp.posthocBF.large.dark.negVSneut@bayesFactor$error * 100,
    N1.peakamp.posthocBF.small.dark.negVSneut@bayesFactor$error * 100,
    N1.peakamp.posthocBF.large.bright.negVSneut@bayesFactor$error * 100,
    N1.peakamp.posthocBF.small.bright.negVSneut@bayesFactor$error * 100
  )
}

# summary
compare.N1.peakamp.posthocBF <- data.frame(
  "comparison" = c("large.dark.negVSneut", "small.dark.negVSneut", "large.bright.negVSneut", "small.bright.negVSneut"),
  "nar" = compare.N1.peakamp.posthocBF[, 1], "nar.p.err" = compare.N1.peakamp.posthocBF.perc.err[, 1],
  "med" = compare.N1.peakamp.posthocBF[, 2], "med.p.err" = compare.N1.peakamp.posthocBF.perc.err[, 2],
  "wid" = compare.N1.peakamp.posthocBF[, 3], "wid.p.err" = compare.N1.peakamp.posthocBF.perc.err[, 3]
)

kable(compare.N1.peakamp.posthocBF)
comparison nar nar.p.err med med.p.err wid wid.p.err
large.dark.negVSneut -1.45681 0.000005 -1.76543 0.000010 -2.09072 0.000012
small.dark.negVSneut -1.39664 0.000005 -1.70120 0.000009 -2.02393 0.000011
large.bright.negVSneut -1.45587 0.000005 -1.76442 0.000010 -2.08967 0.000012
small.bright.negVSneut -1.05302 0.000005 -1.33493 0.000006 -1.64322 0.000007

Follow-up contrasts reveal that the effect of emotion does not explain changes in N1 peak amplitude:

  • large size, high contrast: 0.171113;
  • small size, high contrast: 0.182464;
  • large size, low contrast: 0.171285;
  • small size, low contrast: 0.263177.

Exploratory model selection

omit.from.full.model nar nar.p.err med med.p.err wid wid.p.err
5 emo 1.99368 2.89754 2.31218 2.64091 2.60994 3.10225
2 cont:emo 1.75228 2.87886 2.02680 2.66983 2.35592 3.10442
3 emo:size 1.59916 2.93720 1.93184 2.67864 2.20840 3.10644
1 cont:emo:size 1.43732 2.91708 1.70554 2.67589 2.01439 3.11263
4 cont:size -16.66927 2.97567 -16.80553 2.81724 -16.89753 3.25349
7 size -27.15137 2.94878 -27.29704 2.77292 -27.30376 3.26466
6 cont -38.17203 2.95101 -38.39241 2.73704 -38.46348 3.23530

Exploratory analyses suggest that omitting emotion from the full model improves fitting by 10.096435 times. Removing contrast x emotion, size x emotion, and size x contrast x emotion also improves fitting by 7.589796, 6.902212, and 5.504349 times, respectively.
Conversely, omitting contrast x size or size lowers the explanatory value of the resulting model by 1.98861e+07 and 7.16065e+11 times, respectively. Finally, omitting the factor contrast is maximally detrimental, as it would lower the explanatory value of the resulting model by 4.71643e+16 times.

Peak latency

Peak latency: time (in ms) of the negative value larger than the 3 points (~10 ms at 256 Hz sampling rate) on either side of the peak.
Time window: 150-260 ms post-stimulus onset.
Electrode cluster: TP7 P7 P9 TP8 P8 P10.

size cont emo visualfeats N value sd se ci
large dark neutral large dark 40 187.793 19.4473 3.07488 6.21953
large dark negative large dark 40 190.235 23.3587 3.69333 7.47047
large bright neutral large bright 40 202.149 20.1355 3.18370 6.43963
large bright negative large bright 40 210.742 17.4317 2.75619 5.57492
small dark neutral small dark 40 199.316 21.2240 3.35581 6.78777
small dark negative small dark 40 196.484 17.4944 2.76611 5.59498
small bright neutral small bright 40 214.942 30.5002 4.82250 9.75444
small bright negative small bright 40 213.770 32.1088 5.07684 10.26888

for (k in 1:length(scaling.factor)) {

  ### main effects of size and emotion
  N1.peaklat.sizeplusemo.BF <- lmBF(
    value ~ size + emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peaklat.sizeplusemo.BF, # save model
          file = here::here(paste0("EEG/main/models/N1_peaklat_sizeplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size and emotion
  N1.peaklat.sizebyemo.BF <- lmBF(
    value ~ size * emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peaklat.sizebyemo.BF,
          file = here::here(paste0("EEG/main/models/N1_peaklat_sizebyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of contrast and emotion
  N1.peaklat.contplusemo.BF <- lmBF(
    value ~ cont + emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peaklat.contplusemo.BF,
          file = here::here(paste0("EEG/main/models/N1_peaklat_contplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of contrast and emotion
  N1.peaklat.contbyemo.BF <- lmBF(
    value ~ cont * emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peaklat.contbyemo.BF,
          file = here::here(paste0("EEG/main/models/N1_peaklat_contbyemo_BF_", scaling.factor[k], ".rds")))

  ### main effects of size, contrast, and emotion
  N1.peaklat.sizepluscontplusemo.BF <- lmBF(
    value ~ size + cont + emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peaklat.sizepluscontplusemo.BF,
          file = here::here(paste0("EEG/main/models/N1_peaklat_sizepluscontplusemo_BF_", scaling.factor[k], ".rds")))

  ### interactive effects of size, contrast, and emotion
  N1.peaklat.sizebycontbyemo.BF <- lmBF(
    value ~ size * cont * emo + participant,
    subset(data.peaks, component == "N1" & var == "peak.lat"),
    whichRandom = "participant",
    rscaleFixed = scaling.factor[k],
    rscaleRandom = "nuisance",
    rscaleCont = "medium",
    method = "simple",
    iterations = niter
  )

  saveRDS(N1.peaklat.sizebycontbyemo.BF,
          file = here::here(paste0("EEG/main/models/N1_peaklat_sizebycontbyemo_BF_", scaling.factor[k], ".rds")))

}
compare.N1.peaklat.BF <- matrix(NA, 6, length(scaling.factor))
compare.N1.peaklat.perc.err <- matrix(NA, 6, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {

  ### load models
  N1.peaklat.sizeplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peaklat_sizeplusemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.peaklat.sizebyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peaklat_sizebyemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.peaklat.contplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peaklat_contplusemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.peaklat.contbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peaklat_contbyemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.peaklat.sizepluscontplusemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peaklat_sizepluscontplusemo_BF_", scaling.factor[k], ".rds"))
  )

  N1.peaklat.sizebycontbyemo.BF <- readRDS(
    here::here(paste0("EEG/main/models/N1_peaklat_sizebycontbyemo_BF_", scaling.factor[k], ".rds"))
  )

  compare.N1.peaklat.BF[, k] <- c(
    N1.peaklat.sizeplusemo.BF@bayesFactor$bf,
    N1.peaklat.sizebyemo.BF@bayesFactor$bf,
    N1.peaklat.contplusemo.BF@bayesFactor$bf,
    N1.peaklat.contbyemo.BF@bayesFactor$bf,
    N1.peaklat.sizepluscontplusemo.BF@bayesFactor$bf,
    N1.peaklat.sizebycontbyemo.BF@bayesFactor$bf
  )

  compare.N1.peaklat.perc.err[, k] <- c(
    N1.peaklat.sizeplusemo.BF@bayesFactor$error * 100,
    N1.peaklat.sizebyemo.BF@bayesFactor$error * 100,
    N1.peaklat.contplusemo.BF@bayesFactor$error * 100,
    N1.peaklat.contbyemo.BF@bayesFactor$error * 100,
    N1.peaklat.sizepluscontplusemo.BF@bayesFactor$error * 100,
    N1.peaklat.sizebycontbyemo.BF@bayesFactor$error * 100
  )
}

compare.N1.peaklat <- data.frame(
  "model" = c("size + emo", "size x emo", "contr + emo", "cont x emo", "size + cont + emo", "size x cont x emo"),
  "nar" = compare.N1.peaklat.BF[, 1], "nar.p.err" = compare.N1.peaklat.perc.err[, 1],
  "med" = compare.N1.peaklat.BF[, 2], "med.p.err" = compare.N1.peaklat.perc.err[, 2],
  "wid" = compare.N1.peaklat.BF[, 3], "wid.p.err" = compare.N1.peaklat.perc.err[, 3]
)

compare.N1.peaklat <- compare.N1.peaklat[order(compare.N1.peaklat$med, decreasing = TRUE), ]

kable(compare.N1.peaklat)
model nar nar.p.err med med.p.err wid wid.p.err
5 size + cont + emo 26.6594 1.16392 26.01855 1.24868 25.21809 1.33281
3 contr + emo 24.2288 1.01534 23.82961 1.04759 23.31993 1.08492
4 cont x emo 22.7133 1.26208 21.98327 1.30426 21.12926 1.37474
6 size x cont x emo 21.0690 2.39584 19.23063 2.65336 17.18204 3.00888
1 size + emo 11.1017 1.15394 10.51553 1.20461 9.87124 1.21072
2 size x emo 10.1246 1.40983 9.21805 1.47848 8.28491 1.51659

When using a JZS prior with scaling factor r = 0.707 placed on standardized effect sizes, the model size + cont + emo ought to be preferred.
The best model (size + cont + emo) explains the observed data 8.925722 times better than the second best model (contr + emo).

Paired comparisons

N1.peaklat.posthocBF <- 
  subset(data.peaks, component == "N1" & var == "peak.lat") %>%
  select(participant, condition, value) %>%
  split(.$condition)

for (k in 1:length(scaling.factor)) {

  ### large, dark, negative vs. neutral
  N1.peaklat.posthocBF.large.dark.negVSneut <- ttestBF(
    N1.peaklat.posthocBF$`negativeLargeDark`$value,
    N1.peaklat.posthocBF$`neutralLargeDark`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.peaklat.posthocBF.large.dark.negVSneut, # save model
          file = here::here(paste0("EEG/main/models/N1_peaklat_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### small, dark, negative vs. neutral
  N1.peaklat.posthocBF.small.dark.negVSneut <- ttestBF(
    N1.peaklat.posthocBF$`negativeSmallDark`$value,
    N1.peaklat.posthocBF$`neutralSmallDark`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.peaklat.posthocBF.small.dark.negVSneut,
          file = here::here(paste0("EEG/main/models/N1_peaklat_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds")))

  ### large, bright, negative vs. neutral
  N1.peaklat.posthocBF.large.bright.negVSneut <- ttestBF(
    N1.peaklat.posthocBF$`negativeLargeBright`$value,
    N1.peaklat.posthocBF$`neutralLargeBright`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.peaklat.posthocBF.large.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/N1_peaklat_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds")))

  ### small, bright, negative vs. neutral
  N1.peaklat.posthocBF.small.bright.negVSneut <- ttestBF(
    N1.peaklat.posthocBF$`negativeSmallBright`$value,
    N1.peaklat.posthocBF$`neutralSmallBright`$value,
    mu = 0,
    paired = TRUE,
    method = "simple",
    iterations = niter,
    rscale = scaling.factor[k]
  )

  saveRDS(N1.peaklat.posthocBF.small.bright.negVSneut,
          file = here::here(paste0("EEG/main/models/N1_peaklat_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds")))

}
compare.N1.peaklat.posthocBF <- matrix(NA, 4, length(scaling.factor))
compare.N1.peaklat.posthocBF.perc.err <- matrix(NA, 4, length(scaling.factor))

### model comparison across scaling factors
for (k in 1:length(scaling.factor)) {

  ### load models
  N1.peaklat.posthocBF.large.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_peaklat_posthocBF_large_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  N1.peaklat.posthocBF.small.dark.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_peaklat_posthocBF_small_dark_negVSneut_", scaling.factor[k], ".rds"))
  )

  N1.peaklat.posthocBF.large.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_peaklat_posthocBF_large_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  N1.peaklat.posthocBF.small.bright.negVSneut <- readRDS(
    here::here(paste0("EEG/main/models/N1_peaklat_posthocBF_small_bright_negVSneut_", scaling.factor[k], ".rds"))
  )

  ### model comparison
  compare.N1.peaklat.posthocBF[, k] <- c(
    N1.peaklat.posthocBF.large.dark.negVSneut@bayesFactor$bf,
    N1.peaklat.posthocBF.small.dark.negVSneut@bayesFactor$bf,
    N1.peaklat.posthocBF.large.bright.negVSneut@bayesFactor$bf,
    N1.peaklat.posthocBF.small.bright.negVSneut@bayesFactor$bf
  )

  # percentage of error
  compare.N1.peaklat.posthocBF.perc.err[, k] <- c(
    N1.peaklat.posthocBF.large.dark.negVSneut@bayesFactor$error * 100,
    N1.peaklat.posthocBF.small.dark.negVSneut@bayesFactor$error * 100,
    N1.peaklat.posthocBF.large.bright.negVSneut@bayesFactor$error * 100,
    N1.peaklat.posthocBF.small.bright.negVSneut@bayesFactor$error * 100
  )
}

# summary
compare.N1.peaklat.posthocBF <- data.frame(
  "comparison" = c("large.dark.negVSneut", "small.dark.negVSneut", "large.bright.negVSneut", "small.bright.negVSneut"),
  "nar" = compare.N1.peaklat.posthocBF[, 1], "nar.p.err" = compare.N1.peaklat.posthocBF.perc.err[, 1],
  "med" = compare.N1.peaklat.posthocBF[, 2], "med.p.err" = compare.N1.peaklat.posthocBF.perc.err[, 2],
  "wid" = compare.N1.peaklat.posthocBF[, 3], "wid.p.err" = compare.N1.peaklat.posthocBF.perc.err[, 3]
)

kable(compare.N1.peaklat.posthocBF)
comparison nar nar.p.err med med.p.err wid wid.p.err
large.dark.negVSneut -1.291370 0.000005 -1.588896 0.000008 -1.907163 0.000009
small.dark.negVSneut -1.210616 0.000005 -1.502803 0.000008 -1.817672 0.000008
large.bright.negVSneut 0.503694 0.000002 0.313551 0.000001 0.066594 0.000000
small.bright.negVSneut -1.445018 0.000005 -1.752844 0.000010 -2.077628 0.000012

Follow-up contrasts reveal that the effect of emotion does not explain changes in N1 peak latency:

  • large size, high contrast: 0.204151;
  • small size, high contrast: 0.222506;
  • large size, low contrast: 1.368275;
  • small size, low contrast: 0.173281.

Exploratory model selection

omit.from.full.model nar nar.p.err med med.p.err wid wid.p.err
5 emo 1.844583 3.24770 2.26548 3.67449 2.56387 3.76734
4 cont:size 1.761787 3.32271 2.05090 3.58470 2.43566 3.84750
2 cont:emo 1.467660 3.28640 1.84285 3.57786 2.14626 3.79942
1 cont:emo:size 1.303180 3.25963 1.70487 3.64869 1.99009 3.91498
3 emo:size 0.864005 3.37079 1.21360 3.76316 1.52784 3.96675
7 size -2.447698 3.42704 -2.20890 3.66095 -1.89773 4.08988
6 cont -15.750420 3.56978 -15.66128 3.82710 -15.40421 4.17324

Exploratory analyses suggest that omitting emotion from the full model improves fitting by 9.6358 times. Removing size x contrast, contrast x emotion, size x contrast x emotion, and size x emotion also improves fitting by 7.774863, 6.314525, 5.50065, and 3.365584 times, respectively.
Conversely, omitting size lowers the explanatory value of the resulting model by 9.105681 times. Finally, omitting the factor contrast is maximally detrimental, as it would lower the explanatory value of the resulting model by 6.33296e+06 times.

Interpolated Channels

##   P1.interp.min P1.interp.max N1.interp.min N1.interp.max EPN.interp.min EPN.interp.max LPP.interp.min LPP.interp.max
## 1             0             1             0             2              0              2              0              2

The interpolated channels were mostly identifed outside of the clusters selected for the ERP components defnition (max interpolated channels in clusters: 2); therefore, any potential distortions of the EEG signal due to interpolation was negligible.